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ABSTRACT 


Improvements have been made in the finite element model of the acoustic 
radiated field from a turbofan engine inlet in the presence of a mean flow. The 
problem of acoustic radiation from a turbofan engine inlet is difficult to model 
numerically because of the large domain and high frequencies involved. A numer- 
ical model with conventional finite elements in the near field and wave envelope 
elements in the far field has been constructed. By employing an irrotational mean 
flow assumption, both the mean flow and the acoustic perturbation problem have 
been posed in an axisymmetric formulation in terms of the velocity potential, 
thereby minimizing computer storage and time requirements. The finite element 
mesh has been altered in search of an improved solution. The mean flow prob- 
lem has been reformulated with new boundary conditions to make it theoretically 
rigorous. The sound source at the fan face has been modeled as a combination 
of positive and negative propagating duct eigenfunctions. Therefore, a finite ele- 
ment duct eigenvalue problem has been solved on the fan face and the resulting 
modal matrix has been used to implement a source boundary condition on the fan 
face in the acoustic radiation problem. In the post processing of the solution, the 
acoustic pressure has been evaluated at Gauss points inside the elements and the 
nodal pressure values have been interpolated from them. This has significantly 
improved the results. The effect of the geometric postion of the transition circle 
between conventional finite elements and wave envelope elements has been studied 
and it has been found that the transition can be made nearer to the inlet than 

previously assumed. 
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I. INTRODUCTION 


A. BACKGROUND 

The two predominant sources of noise in a modern jet engine are the jet 
noise and the fan noise. While jet noise provides a challenging opportunity for 
research involving highly non-linear thermal and turbulence effects, the study of 
fan noise has become equally important with the use of high bypass ratio turbofan 
engines in civil aircraft. Although the bypass design has considerably reduced the 
intensity of jet noise by lowering the jet velocity, there is a significant forward sound 
propagation from the fan. This forward radiated acoustic field propagates through 
the inlet duct to be radiated to the far field. Therefore the acoustic analysis of this 
forward propagating noise involves the noise generation by the fan, propagation 
inside the inlet cowling and radiation from the inlet to the far field. This radiated 
acoustic field is highly directional in character with its directivity dependent upon 
the frequency, the mode of the internally propagated acoustic waves, the inlet 
geometry and the mean flow in and around the inlet. The purpose of this work is 
to improve the modeling of the far field acoustical radiation from jet engine inlets 
in low Mach number flows using the finite element method. 

The mathematical modeling of the radiation of turbofan generated noise is 
complicated by the fact that the wavelength of sound radiated may be much 
smaller than the characteristic dimension of the inlet. To resolve the variation in 
acoustic properties near the inlet, it is obvious that a fine mesh must be generated. 
Further complications arise from the fact that sound is radiated to an infinite 
domain and also because of the presence of a mean flow around the inlet. 

The geometries involved in a turbofan engine inlet do not permit a classical 
analytical closed form solution, except when extremely simplified assumptions are 
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made. A wide range of numerical techniques, which have been successfully used 
in duct acoustics, can be thought of as an alternative. These include the Finite 
Element Method (FEM), the Method of Weighted Residuals (MWR), and the 
Finite Difference Method (FDM). Based on previous experience, the FEM showed 
promise in the acoustic radiation problem, provided there was a way to model the 
infinite domain effectively. Fortunately, there has been much work to date in duct 
acoustics which relates to the generation and propagation of noise inside the duct. 
The propagation of sound in the external region, which extends from the throat of 
the inlet to the far field, presents a more challenging computational problem and 
involves the imposition of proper radiation type boundary conditions at a finite 
but distant boundary. In this investigation, the radiation boundary condition is 
imposed by the use of wave envelope elements in the far field. 

B. LITERATURE REVIEW 

The theory relevant to this problem spans many fields in acoustics. Rather 
than discussing each one of them, this section emphasizes the contribution in the 
areas which have a direct link with the problem. Special focus is on duct acoustics 
and the inlet/radiation problem. 

In most practical problems concerning the propagation of sound in ducts, no 
analytical solution is possible, unless extremely simple geometries and assumptions 
are considered. Therefore, the three main numerical techniques which have been 
used in analysis are the Finite Element Method (FEM), the Method of Weighted 
Residuals (MWR) and the Finite Difference Method (FDM). 

Finite element methods have been sucessfully used to model various problems 
in duct acoustics. Astley and Eversman [l], and Eversman, Astley and Thanh [2] 
studied area variation in two dimensions and axisymmetric ducts with different 
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FEM formulations and compared them with MWR results. Majjigi, Sigman and 
Zinn <3 ! used various types of finite elements in the study of simple hard walled 
acoustically treated ducts and compared with results produced by a FDM for- 
mualtion. Various FEM methods have been used to solve the acoustics problem 
with a mean flow. Tag and Lumsdaine [4] used a formulation based on velocity 
potential to save disk storage, while Baumeister [5] indicated that the assumption 
of irrotational acoustic perturbation is valid only for an irrotational mean flow. 

The Method of Weighted Residuals techniques in duct acoustics has been 
useful in situations where appropriate basis functions are obtainable such as hard 
walled or lined ducts with or without flow. Eversman, Astley and Thanh [2] 
compared the results of MWR and FEM methods as mentioned earlier. Eversman 
and Astley (6] investigated the accuracy of MWR compared to exact calculations 
of acoustic transmission based on a one dimensional model for nonuniform ducts 
containing high speed subsonic flow. 

Finite difference methods have not been used extensively in acoustics. How- 
ever, in non-linear problems, for simplified one dimensional models, it seems to 
have a distinct advantage over MWR and FEM. Walkington and Eversman [7] 
studied shocked acoustic waves with a one dimensional model using FDM meth- 
ods. Walkington [8] proposed several schemes to formulate such problems, but 
suggested that extension of the non-linear problem to higher dimensions would be 
difficult. 

The following discussion highlights previous investigations carried out in the 
field of inlet acoustics and the radiation problem. Most work involving inlet acous- 
tics has been experimental in nature, however recently some numerical compar- 
isons have been made. Ville and Silcox [9], and Silcox [10] presented experimental 
results for some standard inlets used by NASA for different flow and geometry 
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configurations. Several analytical and numerical methods have been proposed to 
solve the inlet acoustics problem. Ray acoustics theory (for example Kempton and 
Smith [1 lj) has been combined with numerical flow solutions to analyze various 
inlet geometries. Meyer, Bell and Zinn [12] considered inlet shape and liner de- 
sign by computing far field directivities by an integral method. They numerically 
solved a Helmholtz equation and made an effort not to decouple the far field and 
inlet solutions. Meyer, Daniel and Zinn [13] used the same method as described in 
reference [12] and gave comparisons with experimental results for radiation from 
a pipe and a jet engine inlet. For further improvement, Horowitz, Sigman and 
Zinn used a hybrid FEM-integral technique for cases without mean flow [14], then 
extended this to cases involving mean flow [15], The technique uses a FEM for- 
mulation to analyse the duct interior and then an integral formulation for the far 
field. By guessing a duct exit impedance and solving a duct problem, the far field 
radiation is solved using a Green’s function. The outer boundary impedance is 
then compared to a Sommerfeld condition and the exit impedance is corrected 
iteratively until the results converge. Baumeister [16] used these methods to com- 
pare with experimental data for a JT15D engine under static conditions with a 
low Mach number flow into the inlet. 

The iterative procedure proposed by Zinn et al. turns out to be lengthy and 
costly in terms of both computational tune and storage. To overcome this problem, 
Astley and Eversman [17] employed FEM, wave envelope and infinite element for- 
mulations, and succesfully modeled the sound field for a one dimensional test case 
with no flow. The concept of infinite elements, where the element shape functions 
simulate decay to model an infinite domain, was first proposed by Bettess [18] in 
1977. The application of infinite elements to wave propagation was significant, 
but Astley and Eversman found that “wave envelope” elements, which simulate 
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wavelike behavior in their interpolation functions, model far field acoustic radia- 
tion better than infinite elements. In [19] they laid the foundation for using wave 
envelope elements for inlet radiation in the presence of a flow. They emphasized 
that the use of wave envelope elements relies on the assumption that in the far 
field the sound field approximates that produced by a point source. This allows 
for a coarse mesh in the far field thereby drastically reducing the computational 
time and storage requirements. Astley [20] then validated the concept with simple 
test cases. 

This work is an extension to that done by Eversman, Parrett, Preisser and 
Silcox [21], where they have presented several contributions to finite element mod- 
elling of acoustic radiation from turbofan inlets. This included the use of a tech- 
nique combining finite elements in the near field and wave envelope elements in 
the far field. The use of a frontal solution scheme of Irons [22] resulting in dras- 
tic reduction of in-core storage was also significant. The numerical results were 
verified by comparison with both model scale [10] and full scale [23] test data. 

<? ENHANCEMENTS OF THE TURBOF AN FINITE ELEMENT MODEL 

The finite element model of the jet engine inlet developed in [21] had some 
shortcomings which have been addressed in this study. The finite element mesh 
of the original model had elements whose aspect ratios 1 were not properly main- 
tained. The zone of conventional finite elements outside of the nacelle was orig- 
inally generated in two adjacent regions. This was found to be superfluous and 
was therefore reduced to a single region. In this single region, conventional eight- 
node isoparametric finite elements were generated with their radial thicknesses 

‘Aspect ratio of a two dimensional finite element is the ratio of any two adjacent sides of an 
element. Rule of thumb says that it should not be more than 4:1 to be on the safe side. 
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increasing in geometric progression, thereby maintaining aspect ratios of the ele- 
ments below the safe allowable value. 

The time invariant mean flow problem of the model has a Laplace s differential 
equation with Neumann boundary conditions and the previous solution technique 
to this problem was not theoretically rigorous. The boundary conditions of the 
problem have been modified and a proper solution technique has been incorpo- 
rated. 

A finite element duct eigenvalue problem has been solved on the fan face mesh 
and the resulting finite element modal matrix has been used to model the acoustic 
potential at the fan face boundary as a combination of incident and reflected 
uniform duct eigenfunctions. This boundary condition has been implemented in 
the acoustic radiation problem. 

The problem has been set up with both eight and nine-node quadratic isopara- 
metric elements and results from both the cases have been compared. In the post- 
processing of the solution, the acoustic pressure was observed to be discontinuous 
across inter-element boundaries. This is expected because it can be shown that 
the acoustic potential solution is continuous across inter-element boundaries, but 
its derivative is not. In the original model, the acoustic pressure at a node was 
calculated from the four elements sharing that node and the value of the pressure 
at that node was defined as the average of the these four values. This technique 
gave inadequate results. The improved model evaluates acoustic pressure at gauss 
points inside the element and interpolates the pressure from the interior points to 
the nodes. This has resulted in significantly better results. 

The position of the transition circle, which separates the conventional finite 
element region in the near field and the wave envelope region in the far field, was 
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seen to be creating a significant effect on the final results. It was found out that 
the wave envelope elements were not only capable of modeling the far field but 
also the moderately near field outside the nacelle. Therefore, the transition circle 
could be brought in much closer to the inlet than thought before, and this has 
lead to better results with drastic reduction in the number of degrees of freedom. 
This is probably the most significant result of this study. 
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II. THEORETICAL FOUNDATIONS 


A THE PROBLEM 

The problem posed here is that of a stationary turbofan inlet with internal 
and external mean flow. The results obtained are applicable for an observer fixed 
to the inlet. For an observer on the ground, Doppler shift corrections have to be 
applied to obtain the proper results. 

The turbofan inlet is assumed to be axially symmetric, as is the flow field in 
and around the inlet. The acoustic field generated within the inlet and radiated 
to the far field is generally not axially symmetric, but is conveniently expressed in 
terms of Fourier components in the angular coordinate. Therefore, it is appropriate 
to express the inlet geometry and the entire computational domain in a cylindrical 
coordinate system. 

Figure 1 shows the top half of the symmetric inlet geometry in an z-r plane. 
The surface C n is the nacelle. The nacelle is regarded impervious to both steady 
flow and acoustic perturbations. The surface Cj (fan plane) is the one on which 
the sound source, i.e. the turbofan, is defined (it may or may not have a center- 
body). The acoustic pressure field on the fan face is modelled as a combination of 
incident and reflected (positive and negative) uniform duct eigenfunctions. The 
surface C h (baffle surface) is a boundary of the computational region which for 
a completely accurate representation would be the negative x-axis. However, in 
order to decrease the size of the domain, and also to avoid modeling the rear of the 
engine, C b is chosen at least 90° past the direction of maximum acoustic radiation. 
It may be thought of as a baffle which would admit flow through, yet interfere 
with the forward radiated acoustic field minimally. The boundary Coo is the outer 
boundary of the computational domain in the far field such that at points on Coo 
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the acoustic field can be viewed locally as propagating plane waves. Thus the com- 
putational domain is made finite by imposition of a suitable radiation boundary 
condition at a distant but finite boundary. The boundary should be several 
duct radii from the inlet. The large domain and the fine mesh near the inlet gives 
rise to very large number of degrees of freedom in the finite element model. The 
boundary conditions will be dealt with in more detail in a later section. 

To minimize the problem of huge data storage, an axisymmetric formulation 
based on velocity potential has been adopted, thereby reducing the number of un- 
knowns to one per node. The assumptions of an inviscid fluid and also that of an 
irrotational flowfield has been made (both in steady flow and acoustic perturba- 
tions). The mean flowfield is assumed to be steady, uniform, of low Mach number 
and parallel to the centerline of the engine. These assumptions are quite valid for 
a jet engine during take off or landing approaches. 

B. THE FIELD EQUATIONS 

This sub-section outlines the derivation of the steady mean flow and acoustic 
field equations. It is assumed that the medium is inviscid and non heat conducting 
and that all the processes are isentropic. The field equations have been derived in 
non-dimensional form. Reference density p T and reference speed of sound c T are 
defined at a large distance from the inlet. The reference length is taken as the 
duct radius R. Pressure is non-dimensionalized by p r c’, density by p r , velocities 
and speed of sound by c r , velocity potential by c r R and time by R/c r . 

The governing equations for the problem may be written in nondimensional 
form as 
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Mass conservation : 


|p- + V • ((j’V 1 ) = 0 


( 1 ) 


Momentum conservation : 

£X_ + (V'«V)V' = -^rVp' (2) 

dt P 

Equation of state : 

v = (;) (o ') 1 (») 

where -y is the ratio of specific heats. The non-dimensional speed of sound can be 
written as 

I 



As stated previously, a velocity potential formulation has been proposed to 
reduce the computational time and storage requirements. Since the flow field is 
assumed to be irrotational, the non-dimensional velocity V' can be related to the 
non-dimensional velocity potential $ by 

V' = V$‘ ( 5 ) 

The mass conservation equation (1) in terms of density and velocity potential 
becomes 

^ + V • (p'V$*) = 0 (6) 

dt 

If it is assumed that the reference conditions are taken to be stagnation con- 
ditions, i.e. |V'| = 0, e = 1, d/dt = 0, the momentum equation (2) in terms of 
density and velocity potential may be cast as 

(c‘)* = 1 - b - 1)1^- + j(V*' • V*’)l (7) 
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or, 


t = [1 - (l - l){^- + i(Vf. W)}]'' 1 ’-" 


( 8 ) 


D. LINEARIZED PERTURBED EQUATIONS 

It is assumed that the acoustic quantities consist of small perturbations su- 
perimposed on a steady mean flow field, so that 

= 4>o + <t> ( 9 ) 

P — Pe+P ( 10 ) 

where o subscript denotes the mean flow field variables and the unsubscripted ones 
represent the acoustic perturbation variables. 

Substitution of equations (9) and (10) in equations (7) and (8) and lineariza- 
tion to first order in acoustic perturbations yields 

Pc = [1 - (11) 

+ ( 12 ) 

c* at 

where, c\ — p 2 - * ^ the local speed of sound in the mean flow. 

Similarly, linearization of the mass conservation equation (6) to first order in 
acoustic perturbations yields 

^ + V • {p 0 V4> 0 + Po V4> + pV4> 0 ) = 0 (13) 

at 

Equation (13) is a linear superposition of the mass conservation equation for 
the steady mean flow and the mass conservation equation for the acoustic flow. 
For steady mean flow 

V • (p 0 V<j> 0 ) = 0 (14) 
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and for acoustic perturbation 

— + V . { Po V<t> + pV4> 0 ) = 0 (15) 

dt 

Of the field variables, the physical quantity that can be measured for any 
comparison between theory and experiment is the acoustic pressure. The acous- 
tic pressure is related to the acoustic potential through the linearized isentropic 
equation of state 1 by 

p=-/>.[^ + (V^.V*)] (16) 

Equations (12) and (15) are the basic field equations for the acoustic pertur- 
bation flow. Equation (12) can be used to eliminate p in equation (15) to yield 
a “generalized wave equation” in 4>. Equation (16) can then be used to find the 
acoustic pressure field from the acoustic potential field obtained from a solution 
of (15). The velocity potential field of the steady mean flow needs to be computed 
using equation (14). 

F. THE FTNTTF, ELEMENT MESH AND W AVE ENVELOPE CONCEPT 

The problems for both the mean flow and the acoustic perturbation have been 
solved on the same mesh using a standard Galerkin finite element procedure. This 
sub-section discusses the finite element mesh and introduces the concept of wave 
envelope elements, with the next section discussing the mesh generating scheme 

in detail. 

As shown in Figure 1, the computational region is divided into two major 
regions for conveniently constructing the mesh. The curve C\ marks the border 

2 The linearized isentropic equation of state is 

P= PC o 
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between these two regions. It is important to note that there is a difference in 
the physics involved in generating the mesh in these two regions. In the inner 
region, inside the curve Cj, a conventional finite element mesh based on quadratic 
isoparametric rectangles has been used. The mesh spacing in the general direction 
of noise propagation should be maintained at 4 to 5 elements per wavelength. The 
mesh spacing across the inlet is made fine enough to resolve the transverse modes 
present. In the outer region between C i and Coat a transition occurs between the 
fine mesh and elements which are several wavelengths long. These outer layers of 
elements are called wave envelope elements. 

The major drawback of using a conventional finite element mesh through- 
out the whole domain is apparent when dealing with realistic frequencies. The 
variations of the shape functions of an eight or nine-node isoparametric element 
are quadratic in the local coordinates within each element. Several elements are 
therefore required to accurately represent a single wavelength variation of the so- 
lution in the radial and angular directions. For realistic frequencies, the typical far 
field wavelength of the acoustic field may be several orders of magnitude smaller 
than the overall dimension of the domain. This would demand a very fine mesh 
in the far field and therefore the number of degrees of freedom would become 
prohibitively large. 

To reduce the dimensionality of the problem, wave envelope elements have 
been used in the outer region as an alternative. It has already been assumed that 
Coo is sufficiently far away from the inlet so that the radiated field will behave 
locally as a plane wave propagating outwards from the origin and normal to Coo. 
Therefore, in the outer region, the inlet is assumed to behave as a stationary 
complex source in a uniform flow. Hence the acoustic field in the outer region is 
assumed to be propagating outward with exponential character e • 1, where 
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Tj r is the frequency and %l)(x,r) is the phase, being constant on constant phase 
surfaces. The form of the constant phase surfaces can be visualized by considering 
a simple source in uniform flow (see Figure 2). Note that should be at a 
distance from the inlet sufficient for the exterior flow field to be uniform. The 
constant phase surfaces are found to be 

—Mi + \fx i + fflr* 


V>(x,r) = 


0 2 


(17) 


where (3 7 = 1 - M J . They are circles (of radius R e ) of the form 


(i - MR C ) 2 + r J = R] 


(18) 


The Mach number M here is the Mach number of the exterior flow. 

The shape/basis functions of a typical wave envelope element in the outer 
region are modified from the usual quadratic form to incorporate the complex 
exponential propagation corresponding to a locally outward travelling wave, and 
the reciprocal decay with distance corresponding to a simple source (the velocity 
potential field of a simple source varies as l/r where r is the radial distance from 
the origin). Since the gross features of the harmonic and reciprocal decay solu- 
tion are incorporated into the shape functions, the elements in the outer region 
are required to resolve only the discrepancy between the actual solution and the 
implied harmonic and amplitude variations included in the shape functions. As a 
result, the wave envelope elements can afford to be several wavelengths long and 
the dimensionality of the problem reduces dramatically. 

The modified shape function of node j in a wave envelope element is 

N] = N&e-*^*-**' ( 19 ) 

3 R 

where, 

R = yjz* + /? 2 r* (20) 
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and 


(21) 


— A1 x 4 - R 
^ = — 

as shown earlier. Here N, is the standard shape function corresponding to the j th 
node. The modified shape function of the wave envelope element assumes a value 
of unity at its corresponding node and zero at all other nodes, thereby preserving 
the fundamental property of basis functions. Since the wave envelope elements 
represent the field generated by a simple source in uniform flow, it is expected 
that the elements will be bounded by lines of constant phase and acoustic ray 
paths from the origin as shown in Figure 3. 
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III. MESH GENERATION SCHEME 


This section discusses the generation of the finite element mesh which plays 
a vital role in the formulation of the problem. Both eight-node and nine-node 
quadratic isoparametric elements have been used for the analysis. The mesh gen- 
eration scheme is more or less the same for both the elements, only the connectivity 
and certain other minor parameters need to be altered. Figure 4 shows an eight 
and a nine-node parent element with the local numbering of their nodes. For the 
convenience of constructing the mesh, the entire computational domain has been 
divided into three regions. Figure 5 illustrates the three regions clearly. Region 
I occupies the interior of the nacelle, region II extends from outside the inlet to 
the boundary C\ and region III (the wave envelope region) extends from C\ to the 
outer boundary C 

A. REGION I 

Due to the complex nature of the acoustic field inside the nacelle, a fine mesh 
is generated in order to resolve the variation in acoustic properties. It is separated 
from region II by a circle which we shall call the highlight circle. The highlight 
circle is drawn from the nacelle tip (also known as the highlight) in such a way 
that its center lies at the point of intersection of the z-axis and a line passing 
through the tip of the nacelle at 45° with the z-axis (see Figure 6). 

The inner surface of the nacelle C n extends from the fan face to the tip of 
the nacelle. The centerline of the inlet geometry extends from the intersection of 
the centerbody curve with the z-axis and the intersection of the highlight circle 
with the z-axis. Three-node quadratic line elements lie along the inner surface of 
the nacelle, the centerbody and the centerline. The coordinates of these nodes are 
given as input to generate the mesh in region I. The number of input nodes on 


18 






20 


Figure 6: Geometry of region I 


the inner surface of the nacelle and on the centerbody and centerline is the same, 
to produce a mesh of generally rectangular elements. 

The input data file for the mesh generation program is prepared by a cubic 
spline interpolation routine. The configuration for each of these curves in the inlet 
geometry (i.e. inner surface of the nacelle, centerbody and centerline) is fed into 
the interpolation program in the form of discrete data points. The program then 
fits a smooth curve through these data points to represent that curve by solving a 
tridiagonal system of equations. Convenient nodal points are then chosen on the 
interpolated curve at any desired fraction of the total length. The node points on 
the centerbody and centerline are generated first. These are followed by the node 
points on the inner surface of the nacelle which are at the same fractional distance 
from the fan face (fraction based on the curve length) as their corresponding node 
points on the centerbody and centerline. This has been done to prevent distortion 
in the mesh. 

The fan face has also been divided into several elements not necessarily of 
equal width, each to be represented by a three-node quadratic line element. Figure 
7 illustrates the meshing scheme in this region. The “vertical” element boundaries 
inside the nacelle are formed by arcs of circles. These arcs are drawn through 
corresponding nodal points on the upper and lower boundaries (for example, the 
fifth node on the nacelle inner surface and on the centerbody and centerline, 
counted from the fan face) with the center of the circles lying on the x-axis. Such 
circles are easily constructed as illustrated by Figure 8. (x,, yi) and (x 2 , y 2 ) are 
coordinates of any two corresponding nodal points on the nacelle inner surface and 
the centerbody and centerline respectively. Then the x-coordinate of the center 
of a circle passing through these two points and having its center on the x-axis is 
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Figure 7: Node generation in region I 



the x-axis is given by 


( 22 ) 


_ ( if - A + vj - y?) 

Z ‘ ~ 2(u - x,) 



Figure 8: Geometry of a circular arc in region I 

Now, to preserve the rectangular mesh, each of these circular arcs should 
have the same number of three noded line elements on them and this should 
equal the number of three noded line elements on the fan face. Therefore, each of 
these arcs is divided into the same number of elements with the same fractional 
length (fraction based on the arc length) as on the fan face. Thus, the nodal 
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elements on each of these circular arcs is determined. The nodal coordinates are 
stored in a topology array AD(I, J, K) where I is the global element number, J is 
the local node number and K = 1 assigns the x coordinate, K = 2 the r coordinate, 
respectively, to the array. Proper connectivity relating the local node to its global 
numbering is also generated and stored in a connectivity array AN{I,J), where 
J = element number, J = local node number. This array stores the global node 
number for each node. 

The global node numbering goes from top to bottom of each of the circular 
arcs (see Figure 9) starting from fan face onwards. The element numbering is also 
down each column of elements between adjacent circular arcs and sequenced from 
fan face onwards (see Figure 9). 

B. REGION II 

The mesh in region II becomes polar in nature essentially because of the 
configuration of the domain outside the inlet duct. Region II is separated from the 
wave envelope region III by a constant phase circle, as described previously, whose 
x-intercept is given as input. The constant phase circles are expanding radially 
with the local speed of sound (c) and their centers are moving away along the z- 
axis with the speed of uniform exterior flow (l/) (see Figure 10). This phenomenon 
is very similar to the successive circles of outward ripples created on the surface of 
still water when a pebble is thrown into it. The only difference is that in still water 
the centers of the successive layers of outward moving circular ripples coincide and 
here the centers of the constant phase circles move at a constant velocity. 

From Figure 10, we obtain the equation of a constant phase circle (of radius 
R c ) displaced along the positive x-axis with velocity U (positive direction of U is 
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indicated in Figure 10) 


( 23 ) 


(i - Utf + r 2 = R] 
where R e = ct is the radius of the circle Therefore, 

(x - t/— ) 2 + r 2 = R] (24) 

c 

or, 

(x - MRc ) 2 + r 2 = R\ (25) 

where M is the Mach number of the uniform exterior flow 3 . By setting r = 0 in 
equation (25) we obtain the x-intercept of the circle 

x=(M± 1 )R C (26) 

The positive sign corresponds to the x-intercept on the positive x-axis, 

x = (1 + M)R e 

while the negative x-axis corresponds to the x-intercept on the negative x-axis, 

x = -( l-M)R e 

The circles can be expressed in polar coordinates R and 0 by, 

(Rcos6 — MR C ) l + (iZsmfl) 2 = R\ (27) 

Solving for R in terms of 6 yields, 

R = R e \yJ 1 - M J + (Afcosfl) 1 - McosB] (28) 

Hence, the radial distance R at every angular position 6 on the outer boundary of 
region II is known. 

3 Equation 25 is similar to equation 18 in Section II. 
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The outer surface of the nacelle (portion of the nacelle surface C n starting from 
the baffle C k to the highlight) which forms a part of the inner boundary of region 
II has three-node quadratic line elements along it. The input data preparation for 
this region is the same as in region I. The configuration of the curve representing 
the outer surface of the nacelle is fed into the cubic interpolation routine in the 
form of discrete data points. The program then fits a smooth curve through them. 
Suitable nodal points are then selected at any arbitrary distance along the curve. 
Since the mesh generation in region I precedes that in this region, the coordinates 
of the three-node line elements lying along the highlight circle arc are known. The 
nodal points on the outer surface of the nacelle and on the highlight circle arc 
serve as input for the mesh generation in region II (see Figure 11). 

In this region and also in the subsequent region III, the nodes have been 
generated on and along the acoustic rays from origin. Since the mesh is polar, 
the angular thickness of the elements increases with radial distance because the 
acoustic rays are radial lines diverging from the origin. To maintain proper aspect 
ratio of the elements in this region, the radial thickness of the elements should 
also increase accordingly along acoustic rays moving away from the origin. Now, 
corresponding to each nodal point on the outer nacelle surface and highlight circle 
arc, an acoustic ray is defined and its point of intersection with the outer bounding 
circle of region II is calculated (see Figure 11). Therefore, the radial distance along 
that ray in region II is known. This radial distance is then divided according to 
the number of elements required along the general direction of noise propagation, 
in geometric progression, from the inner boundary to the boundary C%. From 
elementary mathematics, we know that if r It r„ . . . ,r« are n members of a series in 
geometric progression, then the members are related to each other in the following 


27 



u 

V 



28 


Figure 11: Node generation in regions 11 and III 



way, 

r 2 = cr, 
r 3 = cr 2 

r„ = cr n _i 

where c is the common ratio of increment. So, the last member is related to the 
first member by 

r„ = e n - I r, (30) 

Referring to Figure 12 where an acoustic ray intersects with the two boundaries 
of region II, it is obvious that the first and the last members of the geometric 
progression series, i.e. intersection points on the inner boundary and the outer 
bounding circle C\ respectively, are known. Since the number of elements n in the 
radial direction of region II is an input, the common ratio of geometric progression 
is found out using equation (30), 

( outer bounding circle \ n 

z ; — : — — ) 

inner bounding circle J 

Once the common ratio is known, the successive intervals are found out by multi- 
plication with the common ratio as in equations (29). Hence, the nodal points of 
the line elements along that acoustic ray are located. Geometric progression pro- 
vides a gradual increment in the radial thicknesses of elements which is sufficient 
to a maintain proper aspect ratio. 

The nodal coordinate values are stored in rectangular cartesian form in a 
topology array AD(I, J, K) as mentioned before. The connectivity array AN (I, J) 
is also created . As illustrated in Figure 13, the element numbering in this region, 
continues after region I and goes down each column of elements running from the 
baffle surface to the i-axis. The global node numbering also goes down each side 
of the element columns sequenced from the inner boundary consisting of the upper 
surface of the nacelle and highlight circle arc. 
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Figure 12: Geometric progression in region II 



direction 


Figure 13: Element and node numbering in region II 
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C. REGION III 


Region III which consists entirely of wave envelope elements is bounded by 
the transition boundary C j, the outer boundary Coo, the baffle Cj, and the i-axis. 
The wave envelope elements, as discussed before, are large elements bounded by 
acoustic rays and constant phase circles. The string of elements between any two 
successive constant phase circles is referred to as a wave envelope layer. The input 
for mesh generation in this region is the number of wave envelope layers and the 
x-intercept of the constant phase circles bounding each layer. Using equations (26) 
and (28), the inner and outer radii of the constant phase circles bounding each 
such layer is determined. The mean radius of each layer, which is just the average 
of the inner and outer radii, is also calculated. Since the mesh generation in region 
II is complete at this stage, the three-node line elements (note that a three-node 
line element forms a side of an eight or nine-node isoparametric element) on the 
outer bounding circle C\ of region II have been located completely and their global 
numbering is also known. Therefore, corresponding to each nodal point on C i, an 
acoustic ray is defined (see Figure 11) and thereby its points of intersection with 
the inner, mean and outer radii of each wave envelope layer are calculated. The 
rectangular cartesian coordinate values of these intersection points on which the 
nodes lie are stored in the topology array AD(I,J,K). The connectivity array 
AN{I,J) is similarly calculated as in region II. The element and the global node 
numbering follows after region II and is similar to region II. Since the mesh in 
region II is quite fine and that in region III is coarse, care should be taken to make 
a gradual transition in the size of the elements. 

After the mesh is generated, a connectivity check is performed to ensure a 
proper connection between local and global numbering of nodes and uniqueness 
of nodal coordinate values. 
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TV IDENTIFICATION OF THE BOUNDARY ELEMENTS 

Identification of the boundary elements is necessary for proper imposition of 
boundary conditions in the finite element calculations. After the mesh is con- 
structed in the whole domain, the topology (i.e. nodal coordinate) array and the 
nodal connectivity array for the line elements on each of the boundaries of the 
domain are calculated separately. An element identification array N ETY P E is 
set up and different values are assigned to it for different boundary elements for 
identification purposes. 

The setting up of topology and nodal coordinate arrays for each boundary 
surface is accomplished in several subroutines. The nodal connectivity array for a 
boundary is ANL{I,J), where I = element number from 1 to the number of line 
elements along that boundary, and J = the local node number in a quadratic line 
element. This array defines the global node number of the corresponding node on 
that boundary. The topology array is ADL{I , J, K), where I and J are the same 
as above and K = 1 defines the global i-coordinate value of the node, whereas 
K = 2 defines the r-coordinate value of the node. 

An input and output data description for the mesh generation program has 
been described in the appendix. 

F SOME COMMENTS ABOUT THE FINITE ELEMENT MESH 

The acoustic radiation problem is highly mesh dependent but the mean flow 
problem is not very sensitive to the mesh parameters. Since both of these problems 
have been solved on the same finite element mesh, a mesh conforming to the 
acoustic parameters is desired. One of the important factors governing the mesh 
is the number of elements per wavelength which must always be maintained above 
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a minimum value in the main direction of noise propagation to resolve the fine 
variation in acoustic properties. According to the rule of thumb the minimum 
ratio of the number of elements to the wavelength for quadratic elements should 
be 4 or 5. Here a somewhat crude estimate has been made to evaluate that ratio 
along the main direction of sound propagation. 

Since the nondimensional input frequency q r (it is an input to the problem) 
of the sound source on C/ is known, we obtain a ratio of the effective wavelength 
A e (the wavelength of the sound radiated is altered in the presence of mean flow) 

to the reference duct radius R in the following way : 

__ uR _ 2nR 
■ r e A 

since 

A* = (1 - M) A 


therefore, 


2^,1 wx 

T) r = -j— (1 - M) 


^i = — (l-A^ (31) 

R r\ T 

The Mach number M is positive if directed towards the inlet. Now if the number 
of elements per duct radius length is Nr and A is the average width of an element 
within that length, then 



Therefore, using equation (31), the ratio of the number of elements per unit duct 
radius can be expressed as 


Nr 


Vr 

27r(l — M) 


Nx. 


(32) 


where N\,(= K /A) is the number of elements per effective wavelength. For a 
specified number of elements per effective wavelength (for the elements used here 
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N Xt is the goal), equation (32) can be used to determine the number of elements 
per unit of nondimensional length required. This varies as the flow towards the 
inlet varies, and would generally be highest within the nacelle near the fan face. The 
number of elements in the transverse direction within the nacelle or in the angular 
direction in region II is not as critical and is adjusted to maintain the aspect 
ratio of the elements. Another very important mesh parameter affecting the final 
solution is the geometric position of the outer bounding circle C\ (the transition 
circle) of region II. This is discussed in a later section. 

Since the position of the constant phase circles bordering the wave envelope 
layers are user input, care should be maintained to make a gradual transition 
from the small conventional finite elements to the relatively large wave envelope 
elements. For this, the user should be aware of the radial thickness of the last layer 
of conventional elements along C \ . Since the radial thicknesses of the elements in 
region II have been incremented in geometric progression, the radial thickness of 
the last element in region II on the x-axis is 

A r = (c B - c n -')r 0 

where c is the common ratio of increment, n is the number of elements radially 
in region II and r 0 is the x-intercept of Cj. This information is valuable to the 
user for making a smooth transition from region II to III. An example of a finite 
element mesh with 3441 elements (52 elements along the x-direction in region I, 
40 elements radially in region II and 9 wave envelope layers) and 10624 degrees of 
freedom is given in Figures 14, 15, and 16, where Figures 15 and 16 show the mesh 
in regions I and II in detail. The transition circle C\ starts at a nondimensional 
distance of 3.5 from the origin. In region I, a very fine mesh has been generated 
in the x-axis direction due to the complexity of the acoustic field. In the direction 
of the duct radius the mesh has been made gradually more coarse towards the 
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centerbody since the region quite close to the x-axis does not usually fall in the 
main direction of sound propagation. Outside the nacelle, in region II, the mesh 
is coarser than in region I, but still quite dense in the sector bordered by 30° on 
the lower side and 85° on the upper side. This sector usually corresponds to the 
main direction of sound propagation at moderate frequencies (for example, 15.0) 
and low angular mode numbers (for example, 10). The mesh shown as an example 
will be suitable for frequencies upto 15.0. At frequencies higher than 15.0, the 
mesh must be refined in the radial direction to satisfy the number of elements per 
wavelength criterion. The wave envelope elements in the region III allow us to 
have a very coarse mesh in the far field. 
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Figure 14: Finite element mesh inside the nacelle 
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Figure 15: Finite element mesh in regions I and II 






IV. UNIFORM DUCT EIGENVALUE PROBLEM 


For the radiation problem the eigenvectors representing the acoustic modes in 
a uniform duct axe used to implement boundary conditions at the fan face. They 
are conveniently calculated when the mesh is generated. The formulation of the 
eigenvalue /eigenvector problem is discussed in this section. 

A fWF DIMENSIONAL BOUNDARY VALUE PROBLEM 

The fan face Cj is taken to be at a locally uniform part of the inlet. To specify 
the acoustic potential there, a finite element eigenvalue problem has been solved 
in the inlet duct on the fan face. The eigenvalues and eigenvectors obtained from 
the problem are used later to evaluate the boundary condition at the fan face for 
the acoustic radiation problem. The nondimensional acoustic field equation for 
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Figure 17: Duct geometry 


the axisymmetric duct with uniform mean flow shown in Figure 17 is 



+ M— )V- V^ = 0 

ox 


(33) 
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where M is the local Mach number of the uniform flow. It is assumed that the 
acoustic potential in the case of a harmonic acoustic source on the surface C/, 
with time dependence e""', has axial and azimuthal variations of the form 

4> = 4> r {r)e i{r " t ' m9 - k ’ z) (34) 


where m = uR/c t is the local nondimensional frequency (c, is the local speed of 
sound in the duct), k z is the axial wave number and m is an integer representing 
the angular (spinning) mode number 4 . 

Substitution of equation (34) in (33) yields 

- (-)’} - = o (35) 

dr 2 r dr Th Wi r 

It is important to note that the Mach number M in the above equation (35) is the 

local Mach number at the fan face. 


Define 




(36) 


and substitute in equation (35) to obtain Bessel’s equation 

£ + i£ + ( ,4-£),-o 

dr 2 r dr r* 


(37) 


Since the duct has hard walls, the boundary condition prescribed at the wall 
V4> • n = 0 (n is the unit outward normal vector on the duct wall) when inter- 
preted in the one dimensional case yields, 

# Fnr Circular Duct A circular duct corresponds to the case with no centerbody. 

at r = 0, <f)r is finite (38a) 

4 The factor e"* m# accounts for the spinning acoustical modes generated by steady blade loading 
or by the interaction between the rotor and the exit guide vanes. 
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(386) 


at r — R, 



• For Annular Duct An annular duct corresponds to the case with a centerbody. 


dd> T 

at r — r,, — — = 0 

dr 

n d<t>r 

at r = R, — — = 0 

dr 


(39a) 

(396) 


B. FINITE ELEMENT FORMULATION 

A standard Galerkin finite element formulation has been used with three-node 
quadratic Lagrangian elements. A finite element mesh that fits these elements on 
to the fan face has already been dealt with in the previous section. The differen- 
tial equation (35) and the boundary conditions (equations (38) and (39)) which 
compose the boundary value problem for both the circular and the annular duct, 
are approximated by a weak form of the boundary value problem for both the 
circular and annular duct. 


Let ip : n — ♦ R be a smooth function where, the domain fl is [0, r e ] for 
a circular duct and [r^, r 0 ] for an annular one. Multiplication of the differential 
equation (35) with the test function ip and integration over the domain yields 

,d}<Pr . 1 d<pr 


L 01 + rsr + _ ^ w<in - 0 


or 


In + = 0 

Integration of the first term in equation (40) by parts, yields 

d<p r , | f d<P r drp ' 


(40) 




-/ -J L ~r rdr + / (^m - r ^)<p r xprdr = 0 (41) 

| r Jn dr dr J n r z 

where |r indicates that the given term is evaluated at the boundary T. This term 

goes is zero for both the circular and annular duct cases. 
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After discarding the boundary term in equation (41) and rearranging the rest, 


the weak formulation is written as 



d<t> r drjj 

dr dr 


+ —-4> r ^)rdr - 
r i 



(42) 


The following weak problem is posed: find a trial function <f) T : fl — » R 9 
equation (42) holds V smooth ^ : fl — ♦ R. <M r ) and V»( r ) are suitable classes of 
functions whose derivatives are square integrable (from H l space). 


A Galerkin finite element approximation has been used with three-node La- 
grangian quadratic elements. Basis functions Ni,N 2t ...,N n are chosen from an 
n -dimensional subspace of H l . Hence, the test and trial functions can be finitely 
approximated as 

0(r) = c,JV,(r) (43) 

Mr) = dMr) ( 44 ) 

where c,’s and d } ' s are suitable scalars 5 . 


Substitution of equations (43) and (44) in (41) yields a finite element matrix 
formulation of the problem 

</„ {Hrit- + N - N ‘) rdr d ‘ - Kl L N ' N ' rir d - = 0 (45) 


K it 




Kij and M t; are the t, j entries of the stiffness [K] and mass [M] matrices respec- 
tively. 

Non-trivial solutions of equation (45) for the vector d are found if A is an 
eigenvalue of the equation 

(1*1 - *IM|)d = 0 


5 Since <p r is being suitably interpolated between the nodes, the dj ’s here represent nodal values 
of the acoustic potential. 
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or 


{\M\~ X \K\ - [/])d = 0 (46) 

Here A = /c^ is the transverse eigenvalue for the m th radial mode, [/] is the identity 
matrix and d is the vector of nodal values of the acoustic potential. 


The calculations of the global stiffness and the mass matrix are carried out 
at the elemental level and assembly is accomplished using proper connectivity of 
the nodes since, 

m = ew 

\M\ = £|M]- 

where n, is the number of elements in the domain. The element stiffness matrix 
[K}‘ and element mass matrix [M J* are given by 

= f N'N'rdr 
*•*1 

where / n is the integral over the element and Nf is the shape function of the i 
node of an element. 


The matrix eigenvalue problem (equation (46)) is then solved using a QR 
solver. Since the problem is of first order, all of the transverse eigenvalues A are 
real. It is interesting to note that the eigenproblem could have been posed with 
the eigenvalue defined as {k z /rn). The resultant system would have been of second 
order and twice as large as the present one. 

The exact analytical solutions to the differential equation (35) are transcen- 
dental functions called Bessel functions of the first kind J m (* m r) and Bessel func- 
tions of the second kind Y m (K m r) of order m, 

d>r{r) = AJ m {K m r) + BK m (/c m r) (47) 
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Bessel functions are oscillating functions whose amplitudes diminish as /c m r 
increases, and the F m (x m r) become unbounded in the limit as « m r — 0. Therefore 
for the circular duct B = 0 and for the annular one both the constants A and B 
are evaluated by boundary conditions. Application of the boundary conditions on 
equation (47) evaluates the transverse eigenvalues «. m from 

r m (K„r) = 0 (48) 

for a circular duct, and 

JUM + CY^r) = 0 (49) 

for an annular duct for the m tk radial mode. 


C. FORMATION OF MODAL MATRIX 

In classical duct acoustics, it is shown that duct modes can be categorized as 
propagating (cut on) or non-propagating (cut-off). Roughly speaking, lower order 
modes propagate and higher order modes are cut-off. Cut-off modes are those 
which carry no acoustic power and are therefore entirely reactive with energy 
trapped near the source. 


Rearrangement of the terms in equation (36) yields the axial wavenumber k t 
explicitly in terms of the frequency iji and the transverse eigenvalue K m 


kt = m 


- m ± j RT-An^r 


1 -M 1 2 * * * 


(50) 


From equation (50), it is apparent that corresponding to each value of the trans- 

verse eigenvalue x m , there are two distinct values of the axial wavenumber one 

representing a positive 6 (or incident) mode and the other representing a negative 


6 A positive mode is a one which propagates or decays in the positive i-direction. 
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(or reflected) mode. The value of the discriminant in equation (50) determines 
whether the mode is propagating or cut-off. When the discriminant is greater than 
zero, the mode is cut-on. Since the axial variation of the acoustic potential has 
been assumed to be of the form e~ ik ‘* , a positive sign in front of the discriminant 
indicates a positive propagating mode, whereas a negative sign indicates a mode 
propagating in the opposite direction. When the discriminant is less than zero, 
it becomes imaginary and the mode is cut-off. The axial wave number for such a 


mode becomes 


kt = m 


-M ± I - - 1 

1 - M 2 


From equation (34) it can be argued that the amplitude of a wave which is cut-off 
varies along the x-axis as where 0 is 


0 = Vi 


1 - M 2 


Since the amplitude should decay with distance from the source, a negative sign in 
front of the discriminant in equation (51) indicates a positive cut-off mode whereas 
a positive sign indicates a negative cut-off mode. 


The positive and negative duct modes, corresponding to a single transverse 
eigenvalue /c m have the same mode shape (i.e. eigenvector). Since the higher 
order modes are increasingly cut-off, and do not contribute much to the acoustic 
propagation, the first few positive and negative modes have been retained in the 
modal matrix. The modal matrix is an NDOF x (NPOS + NNEG) matrix where 
NPOS is the number of positive modes retained, NNEG is the number of negative 
modes retained and NDOF is the number of degrees of freedom in the system. 
Each column in the modal matrix corresponds to a mode. All the retained positive 
modes have been placed first ordered according to the increasing magnitude of the 
eigenvalue followed by the columns which represent negative modes in increasing 
order. Figures 18 and 19 show the first five acoustic duct modes of an annular 
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duct corresponding to angular mode numbers 10 and 20, respectively. Each mode 
shape is a finite element approximation of a combination of Bessel functions of 
the first and second kind as given by equation (47). Note that for a duct with no 
acoustic lining, the positive and negative propagating mode shapes are the same. 
The transverse eigenvalues and the modal matrix resulting from the calculations 
are used to impose the boundary condition on the fan face in the acoustic radiation 
problem. Details of the imposition of the boundary condition are described in a 
later section. 

The above eigenproblem calculations have been done in the mesh generation 
code because the finite element mesh for the problem is generated along with the 
mesh for the whole domain. An input and output data description for it has been 
elaborated in the appendix. 
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V. TIME INVARIANT MEAN FLOW PROBLEM 

A DERIVATION OF THE FLOW PROBLEMS 

Equation (14) describes the steady mean flow around the inlet of the turbofan 
engine. In the general case the fluid is compressible and equation (14) is nonlinear. 
However, under the assumption of sufficiently low flow hlach number, the flow can 
be approximated as an incompressible one, and equation (14) can be approximated 
by 

V 2 4> 0 = 0 (53) 

which is the Laplace equation. The assumption of incompressibility does not 
impose any extra restrictions on the acoustic perturbation flow equation (15) of 
Section II. D. 

In Figure 1, the curves V in the x-r plane correspond to surfaces in the ax- 
isymmetric space around the inlet. The axis of symmetry r„ is not a physical 
boundary of the domain. In the axisymmetric integral formulation of the prob- 
lem, the boundary integral corresponding to this axis (r = 0) vanishes. Therefore, 
no boundary condition needs to be specified. The far field boundary Too is cir- 
cular (spherical in axisymmetric space). Though it is several duct radii from the 
inlet, the flow effects due to the presence of the inlet cannot be assumed negligi- 
ble. The boundary condition on this curve will be discussed later. The nacelle 
r„ and centerbody T c are impervious to flow. The curve T / represents the fan 
face. The curve F* representing the baffle is a pseudoboundary that does not exist 
physically and corresponds to a porous baffle that admits flow through but affects 
the acoustic field to as small an extent as possible. 

Since the differential equation (53) is linear, it can be split up into three 
different problems, each of which can be solved separately and upon employing 
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the method of superposition, the velocity potential of the actual flow field can be 
determined. The velocity potential for the mean flow is decomposed as follows: 

4> 0 = 4>u + 4 >p 

where fa is the flow field due to the external uniform flow field only (without the 
presence of the inlet) and fa is the flow perturbation due to the presence of the 
inlet only. The boundary condition to be applied at the boundaries r» and Too is 
not clear until and unless we formulate the problem in terms of flow perturbation. 
Our aim is to formulate the entire problem in such a way that the fan face and 
the external flow velocity do not become dependent on the perturbations; rather 
they govern it. 

The perturbation velocity potential 4> p is further decomposed into 

<t>p = <t> 1 + 

where <t>\ is the perturbation potential due to inlet flow alone (fan flow effects 
only) and fa is the perturbation potential due to flow to a blank inlet (effect of 
the presence of the jet engine inlet in the external uniform flow) . Therefore, the 
three flow problems may be posed as 

} Problem I This problem represents the perturbation potential field due 
to inlet flow alone. 

V Vi =0 in fl ( 56a ) 

V fa •n = U f on T f ( 56b ) 

Vfa • n = 0 on T n and T c (56c) 

V<h»n = -^r*n on Too ( 56d ) 

Vfa • n = 0 on T* ( 56e ) 
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where U f is the uniform fan face velocity, A\ is a constant to be determined, n is 
the unit outward normal on the boundary and r is the outward radial vector on the 
outer boundary T „ as shown in Fig.l. It is assumed that on the outer boundary 
To, the effect of the flow field is that of a simple source placed at the origin. Hence, 
the velocity perturbation at the outer boundary is assumed to be radially directed 
inwards and inversely proportional to the square of the radial distance from the 
origin 7 . Therefore the boundary condition (56e) at the baffle boundary T*, which 
is a radial ray, is zero and hence it is impervious to flow perturbations. 

2. Problem II This problem represents the perturbation potential field due 
to a flow to a blank inlet. 

VVs = 0 infl (57a) 

V<t> 2 »n = — V^ u «n onT/ (576) 

V^ 2 • n = — V^ u • n on T n and T e (57c) 

V<£ 2 «n = -^T«n on Too (57d) 

• n = 0 on T* (57e) 

where <f> u is the external uniform flow velocity potential, Ai is a constant to be 
determined and r and n are as mentioned before. Here also the flow at the outer 
boundary is assumed to be that of a simple source placed at the origin. The flow 
perturbation is assumed to be radially outwards and varying as 1 /R 2 , where R 
is the radial distance from the origin. As a result, the baffle boundary (equation 
(57e)) again becomes impervious to flow perturbations. 

3. Problem III The uniform external flow field is generated by 

VV U = 0 in n (58a) 

7 The velocity field of a simple source varies as l/J? 2 . 
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V<£ u • n = U a i • n on To, and Tj 


(586) 


where U a is the external uniform flow velocity. 

Problems I and II are boundary value problems with Neumann boundary 
conditions. Solutions to these problems are non-unique if the value of the unknown 
variable is not specified at one point in the domain Cl. The problems also have 
to satisfy a compatibility criterion which balances the flux of flow across different 
boundaries. This criterion fixes the values of the constants A i and A 2 relative to 
the flow parameters and, therefore, they are not arbitrary. 

A linear superposition of the problems I, II and III gives us the overall bound- 
ary value problem of the mean flow 


V s <£o = 0 inn 

(59a) 

V<f> 0 •n = Uf on T/ 

(596) 

V<t> 0 • n = 0 on r„ and r c 

(59c) 

V<£ 0 • n = — — • n 4- ~^ r • n + U„ i • n on Too 

it 

(59d) 

V<t> 0 • n = U 0 i • n on T* 

(59e) 


The flow perturbation effects of the inlet at the outer boundary are small due 
to the distance of the boundary from the inlet. Also it is to be noted that the 
perturbation boundary condition at the outer boundary Too for problems I and II 
(equations (56 d) and (57<f)) tend to balance each other. Therefore, under these 
conditions the superposed boundary condition (59d) on Too can be written ap- 
proximately as V<t>„ • n « U 0 i • n. The superposition of the elementary solutions 
from problems I, II and III is based on the assumption that the outer boundary 
condition is imposed at a large distance from the inlet. This effectively makes U f 
and U 0 independent of each other. For a given value of U a , any value of U, can be 
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chosen once the elementary solutions are available. Variations in U a requires new 
potential flow solutions to be computed because the mesh depends on U a . 


R FINITE ELEMENT FORMULATION 


] Problem I Since the mean flow field is axisymmetric, there is no variation 
of flow variables in the angular direction. Therefore, the test and trial functions are 
independent of the angular coordinate 0. Let 0 be a real valued smooth function 
defined in the axisymmetric domain 0. Multiplication of the Laplace equation 
with the test function rp and integration over the domain yields 

[ VVitMH = 0 (60) 


By using Green’s theorem, it is determined that 



• n dS 


(61) 


where S denotes surfaces of the axisymmetric volume Cl. Since the problem is 
independent of 6 , the volume integral becomes a surface integral in the x-r plane, 


and the surface integral becomes a line integral, so that 

.d<p i dip . d<pidrp. 


[ /(£*£* + rdxdr = [ .nrJT 

Jr /* dx dx dr dr Jr 


(62) 


Incorporation of the natural boundary conditions into equation (62) results in the 
weak form of the problem 

f !,?h?± + ?h.?±)rdxdr~U,l *rJT-A,( *Lr.nr<ff (63) 
J, J* dx dx dr dr Jr, Jr^ R 

Therefore the following weak problem may be posed: find <p\ : fi ♦ R 2 3 equa- 

tion (63) holds V smooth rp : H — ♦ R*. It is to be noted that <P and are suitable 
classes of functions whose derivatives are square integrable (from the space H'). 


A standard Galerkin finite element approximation has been used for the ma- 
trix formulation. Basis functions N lt N t N n have been chosen from a finite 
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dimensional (of dimension n) subspace of H'. Hence, the test and trial functions 
can be finitely approximated as 


t b — c,Ni{x,r) 


(64) 


<t>i = djNjix^) ( 65 ) 

where c/s and d/s are suitable scalars 8 . Substitution of equations (64) and (65) 
in (63) results in the matrix formulation of the problem 

\[ [ tWdNj- dnW) rdxdT ] d =U/ [ NitdT-Aj N t j- r.nnfr 
| JrJSdx dx dr dr J ' , f]j_ 7r ~ R , 

> 1 ^ F 

Kii ' (66) 

where Kij is the i, j entry of the stiffness matrix [K] and F, is the i tk entry in the 

load vector {.T}. 

2 Problem II In a procedure similar to that of problem I, the weak form of 
problem II is 

[ f(^i2t + ^il^t) r dxdr = -f d>V<t> 0 *nrdT + A 2 [ ^r»nrdT (67) 

J r L [ dx dx dr dr 1 Jv R 

where r = T n U T/ U T e and, <j> and rj) are from HK 

As in problem I, a standard Galerkin finite element approximation has been 
used for the matrix formulation. The matrix formulation of problem II yields 


[//. 


,maN i + maN i)rdxdr ] d 

'r dx dx dr dr 1 l 1 


- j r NiV<t> 0 *nrdT + A 2 N^r^urdT 


( 68 ) 


•Since * is being suitably interpolated between the nodes, d/s here imply nodal values of the 
mean flow potential. 
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where K' tJ is the t, j entry of the stiffness matrix [K'\ and F- is the i th entry in 
the load vector {F'}. 


Details of the stiffness matrix and load vector calculations are dealt with in the 
next sub-section. The constants A x and A 2 of the problems are found by imposing 
the compatibility condition which balances the flux across the boundaries. For 
problem I, it balances the flux across the fan face with the flux across the outer 

boundary T ooj 

u ’L, rdT=A 'L^ r,nTdr 


or. 


Uj_ / r , rdT 
/r.^r.nrdr 


( 69 ) 


For problem II, it balances the flux across the nacelle, centerbody and fan 
face with the flux across Too, i-e., 



• n rdT = A 2 



-^-r • n rdY 
R* 


or, 


f r , V<t> 0 •urdT 

/r. ^r.nrdr 


( 70 ) 


C FINITE ELEMENT CALCULATIONS 

The global stiffness matrices and the global load vectors as defined in equa- 
tions (66) and (68), can be written as the composition of the element stiffness 
matrices and element load vectors respectively. Therefore, for example, in prob- 
lem I, we can write 

m = D*‘] 

in = fun 
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where n t is the number of elements in the domain. The element stiffness matrix 
[K e ] and the element load vector {F e } for the problem are given by 


K: > = 


dN' on; dN[ on; 

dx dx dr dr 


) rdxdr 


F, e = U f [ N'rdT-A l f N'-^rvnrdT 
Jr* Jr* K 


where / n< is the surface integral over the domain of the element, / r « and f r ^ are 
line integrals along element boundaries on the fan face and the outer boundary 
respectively, and N* is the shape function of the t th node of the element. 

Finite element calculations are done based on a parent element with local 
coordinates f and rj as shown in Figure 4. The element shape functions N* corre- 
sponding to each node * in the parent element are standard functions and therefore 
known. 

1. Surface Integrals To perform the finite element calculations on a parent 
element, an element map is constructed. The transformation under which each 
element ft* in the mesh is the image of a fixed parent element under a coordinate 
map T t is constructed as 


nodt$ 

T t : x = £ XiNXs, tf) 

i=1 

(72a) 

node i 

T t :r — J2 TiNttS,*) 

(726) 


1=1 


where nodes is the number of nodes on the parent element, and is 8 or 9 depending 
on whether it is an eight or nine-node element. The element 0 e to which T t maps 
the parent element is completely determined by specifying the x, r coordinates 
(x it rj of all the nodal points of H e . Element shape functions N‘{x,r) are simply 
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obtained from standard parent element shape functions N‘{<;,r)) by 


N'{x,r) = ^(f(x,r),ij(x,r)) 


(73) 


The derivatives of shape functions are obtained by the chain rule of differentiation, 

dNi dNi d$ dNi dr) 

dx df dx * dr\ dx 

(74a) 

dNi dNi d$ dNi dr) 

dr d$ dr ' dr) dr 

(744) 

According to the element map, 


dx dNl 

d < M * 

(75a) 

dx n ^‘ dN e k 
dr) ~ £ 14 dr) 

(754) 

dr n ^* dNl 

* * S r ‘ * 

(75a) 

dr _ dNl 

dr) ~ Tk dr) 

( 7S<0 


By using the above relations (equations (72) through (75)), the element stiff- 
ness expression K* } may be expressed as 


f f[x,r) rdrdx = f ff(f, rj) J(s,ri)dsdr) (76) 

J n, * f)| 

where f n> is the domain integral over the parent element and J (f, rj) is the Jacobian 
of the transformation T t given by 


dx dr dx dr 

d$ dr] dr\ dq 


(77) 


A standard 4x4 Gaussian quadrature rule has been used to evaluate the integral. 
It is important to note that the mean flow calculations have been done both with 
and without the wave envelope elements. In one case, no distinction has been 
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made between the elements in regions I and II and the elements in region III. All 
the regions have isoparametric rectangular finite elements. In another case, shape 
functions for the wave envelope elements in region III, differ from the rest in the 
mesh due to the fact that they simulate the inverse square decay behaviour as 
expected in a field due to a simple source. Mathematically a shape function may 
be expressed as 

N™ = N- 0) ( 78 ) 

where N* is the standard shape function of node * at radius r,. The results in 
both cases were virtually identical for the flow velocity we are concerned with. 
For compatibility with the radiation calculations the wave envelope elements were 
retained. 

2. Boundary Integrals Three-node quadratic line elements lie along the 
boundaries and the generation of their topology and their nodal connectivity have 
already been discussed in the mesh generation scheme. The calculations of the 
line integrals are carried out by integrating along those sides of the parent element 
that are mapped onto the sides F of the actual element fl e along which natural 
boundary conditions are prescribed. For definiteness, it has been assumed that the 
side f = 1 has been mapped onto F. The line integrals have been parametrized 
with respect to r). 

The shape functions used for the line integrals are identical to the standard 
shape functions for the three-node Lagrangian line element. Element maps are 
created as discussed before and the elemental arc length is found by 

dr = Vdx 2 + dr 1 
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or, 


dT = 


(79) 




3 

where j is the jacobian of the transfomation of T) onto the arc length parameter 
in the x-r plane. The dot products, W 0 • n for problem I and r • n for problems 
I and II need to be computed at each node on the relevant boundaries to evaluate 
the line integrals. Note that the constants A\ and A 3 need to be evaluated before 
constructing the load vector. 


n THF, SOL UTION PROCEDURE 

All of the boundary conditions are of the Neumann type and the differen- 
tial equation is the Laplace equation. Hence, there is no unique solution to the 
problems unless a reference value of the mean flow velocity is specified at any 
point in the domain. This does not affect the results because we are interested in 
the derivatives of the potential and not in the absolute values of the mean flow 
potential. By penalization, the potential has been made zero at the intersection 
of the boundaries T„ and Too- This penalization has been made at the elemental 
level. When the stiffness matrix of the element which occupies that node at the 
intersection of T* and T^, is calculated, a very large value (1.0el5) is added to 
the diagonal entry in the matrix corresponding to that boundary node. Hence 
the velocity potential at that node is forced to zero after solution. The penalized 
stiffess matrix for that boundary element looks like the following : 

K\\ 

K ml 

. Km 

where m is the penalized node number (local) and £ (1.0e-15) is the penalty param- 


K im ... K 


In 


+ : ••• K„ 


K nm • * * K n 
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eter. As a check, penalization was carried out at a different point in the domain, 
and the solution was found to differ from the previous solution by an arbitrary 
constant only. 

Since the penalization is carried out at the element level, the stiffess matrix 
and the load vector are never stored in assembled form. As each element stiffness 
matrix and load vector is formed, it is written down onto disk along with its nodal 
connectivity. The frontal solution method of Irons [22] has been used to solve the 
algebraic system of equations [K]{<j>} = {F}. The principles of this technique are 
implied by the Gaussian process of forward elmination and back substitution. The 
frontal process alternates between accumulation of element coefficients (assembly) 
and elimination. Whenever an element is assembled its nodes are kept in active 
storage until their elimination. The active in-core storage at a point of time 
depends only on the “frontwidth” (number of active nodes at that time) which is 
much smaller than the dimension of the assembled matrix. This drastic reduction 
in in-core storage is the most important aspect of this scheme. Details of the 
scheme are, however, not discussed here. 

F.. SUPERPOSITION OF THE SOLUTION FRO M THE THREE PROBLEMS 

After the solutions to problems I and II are obtained, they are added to the 
exterior field velocity potential to obtain the overall mean flow velocity potential of 
the flow field. Solution to the problem III is the uniform flow field whose velocity 
potential is given by 

4> u = U„x + C (80) 

where C is any arbitrary constant. Problems I and II have been solved by penaliz- 
ing the velocity potential at the intersection of Cj and Coo to be zero. Therefore, 
in order to be theoretically rigorous, the uniform flow field velocity potential 4> u 
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should be penalized to zero at that point before superposing the three solutions. 
In this process, the solution vector from the three problems have the same datum 
of reference. Hence, the value of the constant C no longer remains arbritary and 
is calculated as 

C = -U a x p (81) 

where x p is the i-coordinate of the penalized node at the intersection of C* and 

Coo. 

The overall mean flow velocity potential is found out by pointwise addition of 
the solution vectors from the three flow problems 

<t>o = 4>\ + <t>2 + 4>u (82) 

The solutions to the problems I, II and III have been obtained by an input of unit 
velocity at the fan face and in the exterior flow field i.e. Uj — 1 and U 0 — 1. 
Therefore, if the fan face flow Mach number and exterior flow field Mach number 
are Mf and M u respectively, the superposed solution is found by 

4> 0 = + M u [4> a + ^u) (82) 


F. RESULTS AND DISCUSSIONS 

The solution to the mean flow problems I and II are the velocity potential 
values at the nodes in the finite element mesh. Contours of constant velocity 
potential in the field have been plotted in Figures 20 and 21. Figure 20 corresponds 
to inlet flow alone (problem I) with a unit velocity on the fan face and Figure 
21 corresponds to flow into a blank inlet (problem II) with a unit far field flow 
velocity. The contour curves for both the problems are more or less parallel to the 
fan face Cj inside the nacelle and they form concentric circles outside the inlet. 
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The equipotential lines are always orthogonal to the streamlines in a flow field, and 
since the velocity on C „ is radially directed for problems I and II, the constant 
potential curves in the far field are supposed to meet C*, tangentially. From 
Figures 20 and 21 it is obvious that this condition is satisfied upto a certain angular 
distance from the x-axis but at high angles, the contour curves do not quite meet 
Coo tangentially. This is probably because the finite element mesh for the acoustic 
radiation problem is also used for solving the mean flow problem. The mesh in 
region III, corresponding to the wave envelope elements for the acoustic radiation 
problem, is quite coarse for the mean flow problem. This may lead to slight 
numerical inaccuracies in the finite element solution in the far field. However, it is 
not of much concern for the present problem because the perturbation potential in 
the far field is small. Therefore, when the two flow fields are superposed with the 
external uniform mean flow, the effect of the flow perturbation in the superposed 
far field is not noticeable. This is apparent from Figure 22 which shows the 
equipotential lines in the superposed field with a fan face Mach number of 0.5 and 
a uniform far field Mach number of 0.3. It is also to be noted from Figure 22 
that on the far field boundary, the flow is almost fully dominated by the external 
uniform mean flow. Therefore the comment made at the end of Section V.A, that 
the perturbation boundary condition on C* for problems I and II (equations (56d) 
and (57d)) balance each other, is very well satisfied. 
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level potential 
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Figure 20: Contours of mean flow perturbation potential due to inlet flow alone 
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Figure 21: 
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number -0.3, Tan face Maclt number -0.5 


VI. ACOUSTIC RADIATION PROBLEM 


A WEAK FORM AND FINITE ELEMENT FORMUL ATION 

The governing equation for the acoustics problem with flow is given by equa- 
tion (15) in Section II. D 

^ + V . (p 0 Vtf + pV<t> 0 ) = 0 
ot 

where p is given by equation (12) of Section II. D, 

[ 2 * + (v*. • v*)] 

c* at 

c„ being the nondime ns ional speed of sound in mean flow. 

Solutions to the above two equations are desired in the case of a harmonic 
source on the fan face Cj with time and angular dependence given by 
where r\ T is the nondimensional input frequency ( t\ t — u>i2/c r , ui is the input 
frequency) and m is the angular mode number. Since a steady state solution is 
sought, temporal derivatives of acoustic variables are replaced by d/dt = ir) r . To 
formulate a weak problem, let us assume that the trial and test functions are of 

the form 

= *(*,r )««-■-"') (83) 

<Hx,r ,*,() = *(x (84) 

The fact that the test function is taken to be the complex conjugate of the trial 
function is consistent with the definition of the inner product of a complex Hilbert 
space. 

Multiplication of the governing equation (15) with the test function ip and 
integration over the domain yields 

f [ 0 ££ + 0 V • (p„v<£ + pV(t> 0 ) 1 dV = 0 (85) 

Jv dt 
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Integration of the divergence term in equation (85) by parts yields 

f _ V0 • {p 0 V<t> + pVJ>o)} dV = - \ xl>{p 0 V<t> + pV<>o) • n dS (86) 
7v dt Js 

where fy and J s are the volume integral and the surface integral respectively in 

the three dimensional space in and around the inlet. Since the trial and test 

functions have been conveniently chosen, the following weak problem is posed in 

the (x, r) domain fl : find 4>{x,r) : fl — > C 2 9 equation (86) holds V smooth 

V»(x,r) : fl — ♦ C 2 where 

ijjll _ VV> • (po V<t> + pV<t>o) = iVrP^ - Po~^H - pud>,z 

dt T 

— pvd>, r —Po<t>,z t/’.i Po$ ir d>,r (87) 

Here u (= d<f> 0 /dx) and v (= d<j) 0 /dr) are the x and r components of the mean 
flow velocities in the domain fl. ()„ and (), r are the derivatives of variables with 
respect to x and r in the field. 

Substitution of p in the right hand side of equation (87) by equation (12) of 
Section II.D yields 

- Vtl> • [po^4> + pVto) = 

dt 

— \(r)l - C ° r ^ -)^t/ ) + tq r c 0 u(^, x rjj — 4>d>, x ) + ir)rC 0 v(<f),r $ - ) 

c l r 

- [cl - «’)<*, X -(Co - r rl>, r + Uv(*, z r +<t>, r 4> ,* )] (88) 

It is to be noted that the right hand side of equation (88) has products of the 
functions 4> and rj> or their derivatives with respect to x and r. Since the functions 
are complex conjugates of one another, the exponential terms of 4> and as in 
equations (83) and (84) cancel out here. Therefore the right hand side of equation 
(88) is an expression in coordinates x and r only. The volume integral in the left 
hand side of equation (86) reduces to a surface integral over the domain fl in the 
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following way 


f g{x, r,6,t) rdrdxdd = [ f(x,r)rdrdxd6 

Jv Jv 

= 2n [ f{x,r) rdrdx (89) 

J n 

Similarly we can take out a factor 2 tt from the right hand side of equation (86) 
and reduce the surface integral to a line integral. 

J h[x,r,6,t) rdOds = J^p(x,r) rd6ds 

= 27 t J p(x,r)rds (90) 

where ds is the elemental arc length. The constant factor of 2ir therefore cancels 
out from both sides in equation (86). 

The trial and the test functions <f> and respectively are from complex H l 
space. In an analogous procedure to that carried out in Section V.B for the 
mean flow problem, the trial and test functions are finitely approximated using 
the standard Galerkin method. The resulting finite dimensional subspace of H l 
is then constructed. The resulting matrix formulation of the problem yields the 
global stiffness matrix expression of the problem 

Kij = j J h(x, r) rdrdx (90) 

where 

fc(*,r) = * |(l. ! - C ^-)N,N, + N, - N.N,,. ) 

C o r 

+ ir) r c 0 v(N t , r Nj - NiNj„ ) - {cl - u 7 )N,, x N } , x 
~ {cl ~ v^Ni, r Nj„ + uv{N it , Nj„ +Ni, r Nj „ )] (91) 

Ni,N 7 , ... ,N n are basis functions of the finite n-dimensional subspace of H' . 

Conventional eight-node quadratic isoparametric elements have been used to 
model the regions I and II and hence the interpolation functions used here were the 
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same as those for the mean flow problem. The outer region III has been modeled 
with eight-node wave envelope elements. As discussed in Section II. E, the wave 
envelope shape functions are constructed on the assumption that at a sufficiently 
large distance from the inlet, the acoustic field approximates that of a point source 
placed at the origin and are of the form 

where l/R is the nature of decay in the acoustic field due to a point source and 
is the equation of constant phase surfaces (refer to Section II.E). The wave en- 
velope elements are not truly isoparametric in the sense that the element maps 
are created using the standard eight-node quadratic shape functions, but the solu- 
tion is interpolated inside the element using the modified “wave envelope” shape 
functions. 

B. ACOUSTIC BOUNDA RY CONDITIONS 

Equation (86) yields a surface integral of form 

J 0(PoW + P^<t>o) • ndS 

which can be transformed to a line integral over the boundaries in the (x, r) domain 
as has been discussed before. The significance of the combination of terms in 
the integrand is clear because it represents the natural boundary terms which 
are generated by the use of the divergence theorem. The following discussion 
investigates what form the integrand [V>(p 0 V^ + pV<f> e ) • n] assumes on different 
boundaries of the computational domain Q. 

1. Far Field Boundary Ceo On the outer boundary C „ in the far field a 
Sommerfeld radiation condition has been applied. Since there is no reflection in 
the far field this condition assumes that on Coo only an outgoing wave exists. In 
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fact, the modeling of region III with wave envelope elements is complete only when 
this radiation boundary condition at the outer boundary is properly implemented. 


At Coo the acoustic field is assumed to be that of a simple source in uniform 
flow and placed at the origin. The potential for a harmonic acoustic monopole 
(simple source) in a uniform x-direction flow of Mach number M can be written 


as 


where 


4> = 

R 


v/1 - MVo 


4* Pc 


R = + (1 - M 2 )r 2 

1 


0 

k 


(1 - M 2 ) 

Th 

Co 


(-Mi + R) 


(92) 

(93) 

(94) 

(95) 

(96) 


f a is the source strength, q, is the source frequency, p a is the mean flow density 
and k is the local wave number. By taking appropriate derivatives of the acoustic 
potential <t> we get 


d<t> 

- ik 

dx 

[(1 — M 2 ) 

d<t> 

-* - | 

dr 

R 


R } R 2 \ 


(-M + -)- 


R 2 


4> 


(97) 


(98) 


By expanding the right hand side terms of equation (12), we obtain the following 
expression for acoustic density in the field 

^ ± T 

(99) 


_Po 

Co 


96 

ik<t> + M— 
ox 


Substituting p in {p 0 ^<t> + P^4>o) using equation (99), we obtain 
(p c V<£ + pV<t> 0 ) = Po 


— (1 - M 7 ) - ikM<f> 

dx 


_ 96 

X + po—T 

dr 


( 100 ) 
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where x and r are the unit vectors along the x and r axes respectively. Using 
equations (97) and (98), equation (100) can be rearranged to yield 

{p 0 V<t> + pV<t> 0 ) = p 0 -f - (1 ~ R r ~ (^)(*x + rr) (101) 

Therefore, the natural boundary condition at C «, representing the Sommerfeld 
radiation condition for a harmonic acoustic monopole in a uniform flow parallel 
to the i- axis is 

{ Po V4> + pV4> 0 ) . n = ^ -ik- ^ * (<*) (n*x + n r r) (102) 

where n, and n r are the i and r components of the outward unit normal n on the 
outer boundary. 

In the far field, the mesh is constructed on the basis of constant phase circles 
in uniform flow and acoustic rays from the origin. Hence the outer boundary 
as shown in Figure 23 is a constant phase circle. By taking xl> as constant in 
equation (95) we obtain the equation for such a circle to be 


(i - +r 2 = rjj 2 

(103) 

So, a constant phase circle is of radius xl> and having origin at x 
From Figure 23 it is obvious that 

= Mrp and r = 0. 

x - 

n x — cosp = 

(104) 

• r 
n r = sinp = — 

V 

(105) 

thus 


, x 7 ~ Mil>x + r 2 
(n,x + n r r) - 

(106) 


With some algebraic manipulation on the right hand side of equation (106) it can 


be shown that 

n z x + n r r = R (107) 
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Therefore the boundary integral on C «, can be expressed as 

j 0(p o V<£ + pV<t>o) • nrds = Po ~ ' ^ " <t>^ rds ( 108 ) 

As R becomes large (as it is for the outer boundary Coo), the second term in the 
integral becomes negligible with respect to the first and the boundary integral may 
be approximated by 

[ ii>[p 0 V<j> + pV<t>o) •urds « I -P 0 ik<t>xprds (109) 

Jr^, Jr <*> 

In the right hand side of the equation (109), the local wave number k{= Vr/c 0 ) can 

be replaced by rj r because at Coo, «<> = 1. It is important to note that this form 

of the integral assumes a “pc” termination i.e. the wave at Coo behaves locally as 
a plane wave. Therefore, it will be incorrect to compute the boundary integral 
given by equation (109) until and unless the outer boundary Coo is quite far away 
from the inlet. 

In the finite element matrix formulation of the weak problem (86), the bound- 
ary integral (109) yields an n x n matrix A given by 

A„ = J r ip 0 kN,Nj rds (110) 

where n is the number of degrees of freedom (= number of nodes) on the outer 
boundary C^. The outer boundary matrix A is therefore transposed to the left 
hand side of the equation (86) and appropriately added to the stiffness matrix 
given by equation (90) using the nodal connectivity of the outer boundary line 

elements. 

2. Baffle Boundary C» In an attempt to reduce the size of the computational 
domain a baffle surface C* has been modeled at an angle a to the i-axis (see 
Figure 24) which would allow flow through it but affect the radiated acoustic 
field minimally. The baffle C t is swept back far enough from the main direction 
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of sound propagation with the hope that it does not interfere with the forward 
radiated acoustic field. In the region III of wave envlope elements, the natural 
boundary condition given by equation (102) is also valid on Cj since it is valid for 
arbitrary surface shapes. We now investigate the value of the term ( n z x + n r r ) on 

C». 


Referring to Figure 24, we observe that the line representing the baffle C* is 
a straight line passing through the origin and at an angle a to x-axis. Therefore 
its equation is 

y(x,r) = r — 7X = 0 (Hi) 


where 7 = tana is the slope of the line. The outward unit normal n to C* can be 
written as 


n = 


Vy 

|Vy| 


1 

n/TTt* 


(-7X + r) 


( 112 ) 


Thus 


(n,x + n r r) 


>/TTr 


:(-7X + r) = 0 


(113) 


So the boundary integral vanishes on C* in the wave envelope region. Since the 
region II is much smaller than region III, the same approximation has been made 


there while computing the boundary integral on C* with the expectation that any 
errors induced will be localised and will not contribute significantly to the forward 
radiated acoustic field. 


3. The Centerline Since the centerline corresponds to r = 0, the boundary 
integral on the centerline provides no contribution for the axisymmetric formula- 
tion, just as in the mean flow problem. 

4. The Nacelle Surface C w and Centerbody The nacelle surface C n and the 
centerbody are impervious to both steady mean flow and acoustic perturbations. 
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Figure 23: Far field mesh geometry 



Figure 24: Baffle geometry 
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Therefore 


V <j> 0 • n = 0 
V<£ • n = 0 


Hence there is no contribution to the boundary integral in equation (86) from the 
nacelle surface and the centerbody. 


5- The Fan Face The sound source at the fan face C/ has been modeled 
in terms of the duct mode amplitudes. The boundary Cj has been taken to be 
at a locally uniform part of the inlet. The acoustic potential field due to it has 
been expressed as a combination of incident and reflected (positive and negative) 
uniform duct eigenfunctions. The eigenvalues and the eigenvectors from the finite 
element duct eigenvalue problem discussed in Section IV have been used to model 
the natural boundary condition on the fan face C /. The acoustic potential can be 
conveniently written as 

4> = Y. tn e "‘* + " Ie "( r ) + I>« e " kr " le "( r ) ( U4 ) 

n=l r= 1 

where <i> + and <f>~ are incident and reflected duct modal amplitudes, k Zn is the axial 
wavenumber corresponding to positive or negative modes given by equation (50), 
N is the number of modes retained in the expansion and e„(r) is the continuous 
duct eigenfunction corresponding to each retained duct mode. Note that the 
eigenfunctions e„(r) are the same for propagation in the positive and negative 
direction. 


On the fan face boundary the integrand of the boundary integral can be 
expressed by using equation (12) as 


( Po V4> + pV<t>o)**=-Po (* ~ M l) ^ ~ 


(115) 
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where M/ (= V<6 0 /c 0 ) is the axially directed local fan face Mach number. Sub- 
stitution of equation (114) in (115) yields 

(p 0 Vd> + pV«g • n = £ ip 0 f(l - M}) + Vr^] ^e-*- r en (r) 

n— 1 L ° 

+ | >. [(l - Mj) *' + u,^] *;e-‘ ; - , e.(rXU«) 

Tl— 1 

where c„ is the nondimensional local speed of sound in flow at the fan face. 

On the fan face (i = 0, i being measured from the fan face) the acoustic 
velocity potential can be conveniently expressed in terms of duct mode amplitudes 
by substituting i = 0 in equation (114) 


<*>={d(r) ••• e *( r ) e ‘( r ) e "( r ) } 


<t>t ) 


<t>» 




(117) 


Similarly equation (115) can be rewritten as 


(A 


,V<£ + pV<£ 0 )«n = t { Q^ei(r) ... a^e N {r) arej(r) ... cc„e N {r) } 


f 

<t>X 


{ J 

(118) 


where 


a* = Po [(! “ M)) ** + T) r —j- 
So the boundary integral on the fan face C f can be cast as 

, \< 

[ \l)[p 0 V<t> + pV^ 0 ) •nrds = i t/’t 0 B} \ 

Jr f Jr ' { K 


(119) 


K 1 


> rds 


where 


{B} = {e,(r)...e w (r) e,(r) . . . e w (r)} 


( 120 ) 


( 121 ) 
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a. Forced Input to the Acoustics Problem The input to the acoustic radi- 
ation problem has been given in terms of specified values of the incident modal 
amplitudes <t > + , the reflected modal amplitudes <t>~ are obtained as a part of the 
solution. To be precise, the acoustic pressure amplitudes are the inputs and 
by using equation (16) in Section II. D, they are related to the velocity potential 
modal amplitudes by 

P* = -iPoVr (l ~ (122) 

In order to use a forced input on the fan face Cj , generalized coordinates have 
been used on the fan face rather than the nodal values of the velocity potential. 
The generalized coordinates used for this problem are the velocity potential modal 
amplitudes <£* . Equation (117) suggests a convenient transformation from the 
nodal values of the velocity potential on C/ to the generalized coordinates by 

# 1 

Wc, = [M| " (123) 

) 

where {4>}c f is a NF x 1 column vector, NF being the number of nodes on the fan 
face, and [M] is the NF x 2 N transformation matrix ( N is the number of incident 
modes and the number of incident and reflected modes are the same). Since (Af) 
is the matrix of acoustic eigenvectors on Cj, the modal matrix resulting from the 
finite element duct eigenvalue problem serves the purpose. 

This transformation is applied element by element on the fan face Cj. In each 
element on the fan face the boundary nodes (nodes which belong to a fan face 
element and lie on the fan face boundary Cj) are transformed to the generalized 
coordinates but the interior nodes (nodes which belong to a fan face element but 
do not lie on the fan face boundary) remain intact. This tranformation is done by 
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creating an element transformation matrix {M\‘ given by the following relation 


fa 

where {<£}' is the NODE x 1 column vector consisting of nodal values of velocity 
potential, NODE being the number of nodes in a fan face element. The column 
vector on the right hand side is partitioned so that the first 2 N generalized co- 
ordinates are the modal amplitudes and the remaining degrees of freedom fa are 
the velocity potential values of nodes which are interior. Note that the element 
transformation matrix [M\‘ will be different for different fan face elements. 

The effect of this transformation is that in the Galerkin formulation of the 
problem the element stiffness matrix \K e \ y which is basically the integral given by 
equation (90) when calculated over an element, corresponding to elements on the 
boundary C/ are expressed in the generalized coordinates given by the column 
vector in the right hand side of equation (124). The transformed fan face element 
stiffness matrices [if,*] are of the form, 

[K e t ] = [M] tT [K e ]{MY (125) 

Here [K«] is a ( NODE - NB + 2 N) x ( NODE - NB + 2 N) square matrix, 
where NB is the number of fan face nodes belonging to one fan face element. It 
is importrant to note that the first 2 N generalized coordinates consisting of the 
incident and reflected modal amplitudes <£* are common to all the elements along 
the fan face boundary C, but, the remaining degrees of freedom fa corresponding 
to interior nodes of a fan face element are different for different elements. 

b. Finite Element Formulation of the Boundary Integral The boundary in- 
tegral on the fan face C f assumes a very convenient form in the Galerkin method 


w = [*]' 
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where the test or the weighting function is the same as the shape function. Under 
this condition the boundary integral given by equation (120) assumes the form 

K \ 


f 0+ } ( <t>+ \ 

( 0(0*5} | >rds = i jf (B}VB> 

P/ l <t>n J r/ l K J 


rds 


( 126 ) 


Now {B} is the row vector consisting of continuous duct eigenfunctions corre- 
sponding to each retained duct mode. The information that we have regarding 
the duct modes is the finite element modal matrix [M] where each column rep- 
resents a duct mode. But such eigenvectors are discrete. We approximate the 
continuous duct eigenfunctions by the discrete finite element eigenvectors in the 
following way 

{B) = {N} \M\ (127) 

where {N} is the row vector of quadratic basis functions Ni,N 2 , .... Nnf corre- 
sponding to the finite element duct eigenvalue problem. Substitution of equation 
(127) in (126) yields 


ij r {B} T {a*B} | <t>n | rds = i\M\ T {N} T {N} rdr j [a*Af] 


( 1 


1 * 


} ( 128 ) 


The right hand side of equation (128) is evaluated element by element on the fan 
face and each of these element integrals yields a 2 N x 2 N square matrix (C|* given 
by 


[C\ 


' = i\M\ tT N‘} T {N‘}rdr^ [cfM]* 


K 




(129) 


The matrix [C j* calculated in one fan face element is then appended to its corre- 
sponding transformed element stiffness matrix [JT/] given by equation (125). Since 
the first 2N generalized coordinates are the 0*, the matrix [C] e is appended to 
the topmost and leftmost 2 N x 2 N block of \Kf\. Thus the non zero boundary 
conditions in the problem are introduced through the element stiffness matrices 
of the boundary elements. 
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c THE SOLUTION PROCEDU RE 


The frontal scheme of Irons [22], adapted for unsymmetric problems, has been 
used to solve the algebraic system of equations 


m 


4 >~ \ = o 

4 >. 


(130) 


The solution procedure is very similar to that implemented for the mean flow 
problem. The first N generalized coordinates corresponding to <f> + are penalized to 
force in the required input values for the acoustic radiation problem in a manner 
similar to the mean flow case. Such penalization has been carried out at the 

elemental level. 


n RF.SUT.TS AND DI SCUSSIONS 

In this section several example numerical results are presented to demonstrate 
the improvements in the finite element model of the acoustic radiated field of the 
turbofan inlet. The numerical data that can be validated by experiment is that of 
the acoustic pressure in the field. The solution of the acoustic radiation problem 
yields the acoustic velocity potential at the nodes of the finite element mesh. The 
solution is then post processed to yield the acoustic pressure at the nodes by 

p = — Po[*>7r^ + (y4>a • (131) 

which is obtained by using equation (16) in Section II.D. From equation (131) 
it is obvious that the i and r derivatives of both the mean flow and acoustic 
velocity potential need to be evaluated at the finite element nodes before the 
nodal acoustic pressure can be calculated. The differential equations governing 
the mean flow and acoustic radiation problem are of second order and therefore 
the finite element solution space is from H\ This implies that though the solution 
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is continuous across interelement boundaries, the derivatives need not necessarily 
have the same property. Once the velocity potential at the nodes of an element 
are known, the x and r derivatives at any point in the element can be obtained by 


8 -± = t^ 

dx dx 

dr dr 


(132a) 


(1326) 


where n e is the number of nodes in an element. Following the above approach, 
the value of the x or r derivative at a node (which is shared by more than one 
element) might be different when evaluated in the different elements sharing that 
node, since the derivative need not be continuous across interelement boundaries. 
In the original model, a simple average of the nodal derivatives from different 
elements sharing a node was performed to obtain an unique value of the derivative 
at that particular node. From Figure 25, it is clear that at node j, which is shared 
by all the four elements, the derivative with respect to x is obtained for example, 
by 




<t> 1 it +4> i , x +^ S «* J r4>* i 


where 4>\ is the derivative with respect to x evaluated at node j within element t. 
Similar calculations are performed to obtain the derivative with respect to r. 


Test runs have been made at different combinations of source frequencies and 
angular mode numbers for an external uniform flow of Mach number -0.3. Only the 
first radial mode, among the incident ones, is present with a unit modal amplitude. 
Since the acoustic pressure varies over a large range, sound pressure level contours 
have been plotted. 

Figures 26 through 31 show the sound pressure level contour plots with the 
nodal information obtained by the averaging technique. It is observed that the 
contour curves representing the main lobe of radiation in the conventional finite 
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nodes shared by 
more than one 
element 



Figure 25: A node shared by all the elements in a 4 x 4 mesh 


element region II are extremely jagged in nature and do not represent the true 
nature of the acoustic field there. But strangely enough, they smoothen out totally 
in the wave envelope region HI. This is true for low frequencies like 12.0, which is 
on the lower side, and also for 20.0 which is on the higher side. Also it is noticed 
that this erroneous behavior is more spread out in the field at higher frequencies. 
Besides this, spurious reflections from the baffle are also present in the radiated 
field especially at high frequencies. 

To improve on the results, nine noded quadratic isoparametric elements were 
used for the analysis instead of the eight noded ones. The presence of the extra 
ninth node at the element center did not provide better results. In fact, the results 
when compared with the eight noded ones were almost identical. Therefore, to 
reduce the dimensionality 9 of the problem, the nine noded elements were discarded. 
However the analysis is almost the same with both kinds of elements, except the 

9 The number of degrees of freedom associated with N eight noded elements in the mesh is less 
than A nine noded elements by N . 
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Figure 27: Sound pressure level contours inside the nacelle and in the near Geld 
from nodal data; tj r = 12.0, m = 10 


vel pressui 
C 100.0 
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level pressure (dB) 
C 135.00 



Figure 30: Sound pressure level contours in the whole domain from nodal data; 



vel pressure (< 
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the nodal connectivity and the coordinate array is generated in a different way 
while constructing the finite element mesh. 
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Figure 32: 2 x 2 grid of Gauss points in the parent element 

For one dimensional line elements, it can be mathematically proved that the 
optimum points in the element for the derivative to be calculated are the ones 
corresponding to the standard Gaussian quadrature points in the parent element. 
Since these points are internal to an element, the derivatives of the solution cal- 
culated there are unique and continuous. Extending this feature to the two di- 
mensional case, we seek to develop a scheme for calculating the acoustic pressure 
at the interior points of the elements corresponding to a standard 2x2 grid of 
Gauss points in the parent element (see Figure 32). 

The nodal velocity potential solutions (both mean flow and acoustic) are in- 
terpolated at the Gauss points in an element using 

>y*) = H <t>oiN;{x„y , f ) (133a) 

i=l 
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( 1336 ) 


n< 

4>{x a ,y a ) = £d>.N'( x a ,y g ) 

t=l 

where (x a ,y a ) are the coordinates of the map of the Gauss points' 0 in the element. 
The derivatives of the potentials are then taken at the Gauss points by appro- 
priately taking the derivatives of the element shape functions at the gauss points 


since 


dd»(j g ,y g ) 




dA7(x„y,) 


(134al 


d<t>{x a ,y a ) P* . dNt^Va) 
-L* dr 


(1346) 


dr 


i=i 


Once the velocity potentials (both mean flow and acoustic perturbation) and their 
derivatives have been calculated at the Gauss point grid inside the element, the 
acoustic pressure is evaluated there using equation (130). The sound pressure level 
contour plots from the data at the Gauss points give smoother curves and represent 
the acoustic radiated field much better. This is reflected in Figures 33 through 35 
where the sound pressure level curves have been plotted for the same combinations 


of source frequencies and angular mode numbers as in Figures 26 through 31, 
but using the data at the Gauss points. Therefore the averaging technique for 
evaluating the nodal pressures as mentioned before seems to be quite inadequate 
for the acoustic radiation problem especially for higher frequencies. The proper 
and theoretically more rigorous way to obtain the nodal acoustic pressures is to 
obtain the pressures at the Gauss points and then interpolate the nodal values 

from them. 


The next part of the study was aimed at investigating the effect of the geo- 
metric position of the transition circle C, in the mesh. If the transition circle is 
very far away from the inlet, then, at realistic frequencies which are usually the 

10 Gauss points are actually in the parent element. They are mapped onto corresponding point, 
in the actual element. 
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higher ones, the number of conventional finite elements radially in region II has 
to be quite large to satisfy the number of elements per wavelength requirement. 

If the transition circle can be brought closer to the inlet the number of elements 
radially in region II needed to satisfy the number of elements per wavelength re- 
quirement will lessen. In such a case, the wave envelope region III becomes larger, 
but the number of wave envelope elements in that region does not have to increase 
proportionally because these elements have inverse decay and exponential terms m 
them to model the field due to a simple source. As the transition circle is brought 
closer to the inlet, the number of degrees of freedom associated with the problem 
drastically reduces, even at higher frequencies. Theoretically it cannot be brought 
very close to the inlet because the field there does not behave as one due to a 
harmonic acoustic monopole in uniform flow. Figures 26 through 35 show sound 
pressure level curves with the transition circle C, at 3.5 duct radius from the inlet. 
Figures 36 and 37 show sound pressure level contours at r) r = 20.0 and m = 20 
with the transition circle at 2.5 duct radius and 1.5 duct radius, respectively, from 
the inlet. Comparing the results in Figures 35, 36 and 37 for r} T = 20.0 and m = 20, 
it is very surprising and also encouraging to observe that the the sound pressure 
level contours get smoother as the transition circle is brought closer to the inlet, 
even as close as 1.5 duct radius from the origin (the i-intercept of the transition 
circle is at a nondimensional distance of 1.5 from the origin). Furthermore, a 
drastic reduction in the number of degrees of freedom occurs. Figure 35 shows 
level curves obtained from a mesh of 4666 elements (14349 dof) while Figures 36 
and 37 show the same level curves obtained from meshes having 3441 elements 
(10624 dof) and 2461 elements (7644 dof), respectively. As the transition circle is 
brought in closer to the inlet the level curves become significantly smoothen along 
with the tremendous reduction in the dimensionality of the problem. Even though 
the transition circle is close to the inlet, the curves have a smooth reflection free 
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transition from the conventional finite element region II to region III. We can infer 
from this that though the wave envelope elements are meant for modeling the far 
field, their shape functions which have the conventional finite element shape func- 
tion expression in addition to the inverse decay and exponential terms, are also 
capable of modeling the moderately near field outside the inlet. As the transition 
circle is brought in closer, the angular resolution of the conventional finite element 
mesh in region II increases. This probably contributes partly to the improvement 
in results. However, in an attempt to reduce the dimensionality of the problem, it 
will be incorrect to apply the u pc n termination (the Sommerfeld radiation bound- 
ary condition) at a boundary quite close to the inlet. Further investigation needs 
to be carried out to understand the phenomenon more clearly. 

Referring to Section VI.B we notice that the acoustic boundary condition at 
the portion of the baffle C* belonging to region II is not properly applied with the 
hope that the errors due to it will be localized. But we observe that especially at 
higher frequencies and higher angular mode numbers a significant level of sound 
is radiated around the lip of the nacelle towards the baffle. The portion of the 
baffle in region II does not act as a reflection free boundary and therefore creates 
incorrect standing wave patterns near the baffle. Spurious reflections from the 
baffle in region II are also observed at lower frequencies, but the intensity is much 
less. As the transition circle C, is brought closer to the inlet, the part of the 
baffle belonging to region II becomes less and therefore the spurious reflections 
reduce significantly in intensity because the part of the baffle belonging to region 
III is reflection free. The reduction in dimensionality of the problem as well as 
improvement in results as the transition circle is brought closer to the inlet are 
significant results of this study. 
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VII. CONCLUSIONS 


Several improvements to the finite element modeling of acoustic radiation 
from turbofan engine inlets have been presented. They are enumerated below: 

. The finite element mesh has been improved in search of an improved solution. In 
particular, the aspect ratios of the conventional finite elements in the moderately 
near field (region II) have been maintained, especially in the main direction of 
sound propagation. The number of conventional finite element regions outside 
the nacelle was reduced to one thereby eliminating superfluous coding and also 
obtaining some reduction in the total number of degrees of freedom in the mesh. 

• The time invariant mean flow problem has been reformulated with new boundary 
conditions and a proper solution technique has been incorporated. 

• A finite element duct eigenvalue problem has been solved on the fan face mesh 
and the resulting modal matrix and the eigenvalues have been used to incorporate 
a source boundary condition on the fan face in the acoustic radiation problem. 

• The acoustic velocity potential at the sound source has been modeled as a 
combination of the positive and negative propagating duct modes evaluated by 
the finite element duct eigenvalue problem. By employing this, a finite element 
formulation of the boundary integral on the fan face has been obtained. 

• In the post processing of the solution, the acoustic pressure was observed to 
be discontinuous across inter-element boundaries. The technique of averaging the 
pressures at a node, calculated from each element sharing that node, was found to 
give poor results especially at higher frequencies. An improved way of evaluating 
the acoustic pressure at the Gaussian quadrature points inside the elements and 
then interpolating it to the nodes has significantly improved the results. 
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• The geometric position of the transition circle bordering the conventional finite 
element region and the wave envelope region was seen to be creating a significant 
effect on the results. It was found out that the wave envelope elements were not 
only capable of modeling the far field but also the moderately near field outside the 
nacelle. Therefore, the transition circle could be brought in much closer to the inlet 
than thought before, and this has lead to better results with a drastic reduction 
in the number of degrees of freedom. This is probably the most significant result 
of this study. 

These contributions have been implemented in computer programs which are 
capable of predicting the radiation pattern at frequencies which are comparable to 
those of actual turbofan engines. With the current version of the computer pro- 
grams it is possible to accomodate 21000 degrees of freedom which would predict 
the radiation pattern fairly well at quite high frequencies (for example, t) r = 35) 
provided the transition boundary circle is not very far away from the inlet. Future 
work should aim at imposing a proper and accurate boundary condition on the 
baffle in the conventional finite element region II in order to make it reflection free. 
The phenomenon of improvement in the results as the transition circle is brought 
closer to the inlet should be investigated in more detail. The representation of 
the sound source in terms of the duct mode amplitudes may not be the best way 
to model it. Alternative ways should be investigated and the results should be 
compared with experimental data to see which model works best. 
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APPENDIX A 


TTSF.F’S MANUAL 

/t DATA DESCRIPTI O N OF THE MESH GENERATION CODE 

1, Input data The finite element mesh generation code PRATMESH not only 
generates the finite element mesh for the problem but also solves the duct eigen- 
value problem on the fan face mesh. Therefore the input data descriptions for these 
two steps are included together. The program is dimensioned for working with 
a maximum of 5000 eight node quadratic isoparametric elements, 21000 nodes, 
150 boundary line elements on the combination of upper and lower surface of the 
nacelle, 50 boundary line elements on the fan face and 125 boundary line elements 
on the far field boundary. The input data file structure has distinct blocks of input 
data referred to as cards. Each card begins on a new line and the input data on 
the card is formatted. The record length is of a maximum of 80 characters. If the 
data format requires more than 80 characters a card is continued on additional 
lines. The input variables should be input in the file in the sequence given. 


Card 

Variable 

Format 

Description 

1 

NLINU 

15 

Number of three-node line elements to de- 
scribe the upper surface of the nacelle 


NLINI 

15 

Number of three-node line elements to de- 
scribe the lower surface of the nacelle 


NY 

i 

15 

Number of elements along the duct radius 
in region I, which is equal to the number of 
three-node line elements on the fan face. The 
fan face is the plane x = constant, at which 
the input duct modal amplitudes are speci- 
fied. 

,i 
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Card 

Variable 

Format 

Description 

1 

NX2 

15 

Number of elements radially in region II 


NWEEL 

15 

Number of wave envelope layers 


NPRINT 

15 

= 0, do not write nodal coordinate array 
= 1, write nodal coordinate array 


NDONT 

15 

= 0, do not write output data file 
= 1, write output data file 


NCNTR 

15 

= 0, there is no centerbody 
= 1, there is a centerbody 


NCNTR1 

15 

Sequence number of the first line element on 
the centerbody. It is the first element from 
the intersection of the centerbody and the 
fan face. 


NCNTR2 

i 

15 

Sequence number of the last line element on 
the centerbody. It is the element at the in- 
tersection of the centerbody and the x-axis. 

2 

PCNT(I) 

6F10.0 

1=1, NY; end node locations of the three- 
node line elements lying along the fan face. 
The node locations are given as fractions of 
the fan face width and starting from the in- 
tersection of the fan face and the lower sur- 
face of the nacelle. The first node therefore 
has a zero fractional distance and is not an 
input. If there are 5 line elements along the 
fan face then a typical input for this array 
would be 0.2, 0.4, 0.65, 0.88, 1.00 
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Variable 


Description 


ADSHPU(I, J, K) 6F10.0 


ADSHPI(I, J, K) 6F10.0 


CBSHPI(I, J, K) 6F10.0 


I = 1, NLINU; J = 1, 3; K =1, 2 
nodal coordinates of the three-node 
line elements defining the upper sur- 
face of the nacelle. I is the ele- 
ment number, J is the local node 
number, K = 1 defines x coordinate 
value and K = 2 defines r coordinate 
value. The elements are sequenced 
from the baffle surface to the to the 
tip of the nacelle. Each card has the 
nodal coordinate information of one 
line element. Therefore, the number 
of records for this card will be the 
number of line elements along the 
upper surface of the nacelle 

I = 1, NLINI; J = 1, 3; K =1, 2 

nodal coordinates of the three-node 
line elements defining the lower sur- 
face of the nacelle. I is the ele- 
ment number, J is the local node 
number, K = 1 defines x coordinate 
value and K — - 2 defines r coordinate 
value. The elements are sequenced 
from the fan face to the tip of the 
nacelle. The number of cards will 
be the number of line elements along 
the lower surface of the nacelle 

I = 1, NLINI; J = 1, 3; K =1, 2 
nodal coordinates of the three-node 
line elements defining the center- 
body and the centerline of the com- 
putational domain. The array is 
similar to that in cards 3 and 4 and 
so also is the format.The elements 
are sequenced from the fan face to 
the intersection of the highlight cir- 
cle with the x-axis. 
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Card 

Variable 

Format 

Description 

6 

R2 

FlO.O 

i-intercept of the outer bounding circle C\ of 
region II in multiples of the inlet duct radius 

7 

RLAYER(I) 

6F10.0 

I = 1, NWEEL 

x-intercepts of the outer bounding circles of 
the wave envelope layers in multiples of the 
inlet duct radius 

8 

MT 

15 

Angular mode number 

NPOS 

15 

Number of positive modes retained in the 
modal matrix 

VMACH 

F10.5 

Freestream flow Mach number outside the 
nacelle (positive directed towards the inlet) 


2. Output data PRAT MESH has two output data files - unit 6 and 20. Data 
file 6 is the printed output data file by default and is well documented in itself. 
Therefore, it is not described here.The output data file 20 contains all the infor- 
mation about the finite element mesh and serves as an input to the finite element 
calculations in subsequent codes. The output data file structure, like input file unit 
5, has records which are classified under the heading of different cards because of 
the varied nature of input parameters. Each card begins on a new line and the 
input data is formatted. The record length is of a maximum of 80 characters. 
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Description 





110 Number of elements in the domain 


110 Number of three-node line elements to de- 
scribe the surface of the nacelle 


HO Number of elements in region I along the 
duct radius 


110 Number of elements in region II or region III 
in the angular direction 


NNXl 110 Number of nodes along z-axis in region I 


NNX2 110 Number of nodes radially in region II 




NNTH1 


NNTH2 




Number of nodes in region I along the duct 
radius 

Number of nodes in region II in the angular 
direction 


NNODE 110 Number of nodes in the domain 


NWEEL 110 Number of wave envelope layers 
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Variable 


Description 

NXl 

no 

Number of elements along z-axis in region I 

NX2 

no 

Number of elements in region II along the 
radial direction 

NY 

no 

Number of three-node line elements along 
the fan face 

NY2A 

no 

Number of three-node line elements along 
the upper surface of the nacelle 

NY2 

no 

Number of elements in region II in the angu- 
lar direction 

NY3 

no 

Number of elements in region III in the an- 
gular direction 

NCNTR 

no 

1 or 0 value deciding the presence or absence 
of the centerbody 

NCNTRl 

no 

Sequence number of the first centerbody el- 
ement 

NCNTR2 

no 

Sequence number of the last centerbody ele- 
ment 

NLC 

no 

Number of centerbody elements 
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Description 


D15.8 x intercept of the outer bounding circle 
Ci of region II 


6 RLAYER(I) 4D15.8 x-intercept of the outer bounding circles 

of the wave envelope layers 




AN (I, J) 7H0 Nodal connectivity array for the elements 

I = element number, J = local node num- 
ber 


8 NETYPE(I) 


ANL2(I, J) 




12 ANL3(I, J) 


I = 1, NE 

identification number of elements 

nodal coordinate/topology array of the el- 
ements 

I = Element number, J = local node num- 
ber, K = x or r coordinte specifier 


Nodal connectivity array for the bound- 
ary line elements along the nacelle surface 
I = line element number, J = local node 
number 


Nodal connectivity array for the bound- 
ary line elements along the fan face 
I = line element number, J = local node 
number 

Nodal connectivity array for the bound- 
ary line elements along the outer bound- 
ing circle Coo of the domain 
I = line element number, J = local node 
number 
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Card 

Variable 

Format 

Description 

13 

ANLC(I, J) 

7110 

Nodal connectivity array for the bound- 
ary line elements along the centerbody 
I = line element number, J = local node 
number 

14 

ADL1(I, J, K) 

4D15.8 

Nodal coordinate array of the boundary 
elements along the nacelle surface 
I = element number, J = local node num- 
ber, K = x or r coordinate specifier 

15 

ADL2(I, J, K) 

4D15.8 

Nodal coordinate array of the boundary 
elements along fan face 
I = element number, J = local node num- 
ber, K = x or r coordinate specifier 

16 

ADL3(I, J, K) 

4D15.8 

Nodal coordinate array of the boundary 
elements along 

I = element number, J = local node num- 
ber, K = x or r coordinate specifier 

17 

ADLC(I, J, K) 

4D15.8 

Nodal coordinate array of the boundary 
elements along centerbody 
I = element number, J = local node num- 
ber, K = x or r coordinate specifier 

18 

PCNT(I) 

4D15.8 

1=1, NNTHl; fractions in which the fan 
face has been divided for input nodes (see 
input data description) 

19 

MT 

no 

Angular mode number 
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Card 


Format 

Description 

19 

NPOS 

110 

Number of positive modes retained in the 
modal matrix 


NNEG 

110 

Number of negative modes retained in the 
modal matrix 

NVECT 

no 

Total number of positive and negative 
modes retained in the modal matrix 

20 

VMACH 

D15.8 

Freestream flow Mach number outside the 
nacelle (positive directed towards the inlet) 

21 

DD(I, J) 

4D15.8 

I = 1, NNTH1; J = 1, NVECT 
Truncated modal matrix from the finite el- 
ement duct eigenvalue problem 

22 

L 

VKAP(I) 

4D15.8 

I = 1, NPOS 

Transverse eigenvalues of the annular duct 


ft data descripti o n op the mf,an FLOW CQD £ 

1 . Input data The time invariant mean flow problem on the finite element mesh 
is solved in the PRATFLOW code. The mesh information in the output file unit 
20 from the PRAT MESH program serves as an input to the PRATFLOW code. 
The user input for this program is the data file unit 5. It has only one card and 
the input data is formatted. The record is a maximum of 80 characters. The mput 
variables should be input in the file 5 in the sequence as described. 
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Card 


Variable 


Format 


Description 


1 NTYPE 15 = 1, calculates the mean flow velocity poten- 

tial only for flow from infinity into the blank 
inlet 

= 2, calculates the velocity potential for inlet 
flow alone 

= 3, calculates both cases sequentially and 
writes the solution vector to disk for use by 
the superposition program 

PRINTl 15 # 0, beginning row and column of the as- 

sembled stiffness matrix to be printed 
= 0, do not print 

PRINT2 15 # 0, final row and column of the assembled 

stiffness matrix to be printed 
= 0, do not print 

NPLOT 15 = 0, contour plotting routine bypassed ma- 

trix to be printed 

= 1, plot contour level curves for the solution 
vector 

NPRINT 15 = 0, do not print nodal coordinate array 

= 1, print nodal coordinate array 


2. Output data The output data files from the PRATFLOW code are files - 
unit 6 and 21. File 6 is well documented in itself and is not described here. File 21 
contains the nodal mean flow velocity potential of problems II and I in that order. 
It serves as an input to the PRATV EL program. In addition to the saved files 6 
and 21, seven unformatted scratch files - units 1, 2, 4, 8, 15, 16 and 17 are used. 
Files 15, 16, 17 are direct access and files 1, 2, 4, 8 are sequential access. Record 
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lengths for the direct access files 15, 16, and 17 are 36, 144 and 36 respectively. 
The number of records in each direct access file is the number of elements in the 

domain. 

q DATA DESCRIP T ION OF THE VELOCITY POTENTIA L 
SUPERPOSIT ION CODE 

1. Input data Problems I, II and III of the time invariant mean flow problem are 
superposed in the PRATVEL program to obtain the mean flow velocity potential 
at the nodes. Subsequent calculations to obtain the nodal mean flow velocity is 
also carried out in this code. The finite element mesh information in the output 
file unit 20 from the PRATMESH program serves as an input to the PRATVEL 
code. The output data file unit 21 from the PRATFLOW code is another source 
of input. It contains the nodal values of the velocity potential of problems II and 
I of mean flow in that order. The user input for the PRATVEL code is data file 
unit 5 which has only one card and the data is formatted. The input variables 
should be input in file 5 in the sequence described. 


Card 

Variable 

Format 

Description 

1 

VMIN 

F10.0 

Average compressible inlet Mach number at 
fan face (positive directed towards inlet) 
based on local speed of sound 


CFS 

F10.0 

Free stream speed of sound, outside the na- 
celle 


RHOFS 

F10.0 

Free stream density of air, outside the nacelle 
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2. Output data The output data files for the superposition program are files unit 
6 and 22. Data file 6 is well documented in itself and hence is not described here. 
File 22 serves as an input data file for the acoustic radiation program. It contains 
some flow parameter values and the nodal values of the mean flow velocity. The 
output data description in file 22 is given below. Each card of data begins on a 
new line. 


Card 

Variable 

Format 

Description 

1 

VMIN 

D15.8 

Average compressible inlet Mach number at 
fan face (positive directed towards inlet) 
based on local speed of sound 


CFS 

D15.8 

Free stream speed of sound outside the na- 
celle 


RHOFS 

D15.8 

Free stream density of air outside the nacelle 


CSTAG 

D15.8 

Stagnation speed of sound 

2 

RHOSTG 

D15.8 

Stagnation density 


TZERO 

D15.8 

Stagnation temperature 


CF 

D15.8 

Speed of sound at fan face 


RHOF 

D15.8 

Density at fan face 
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Card 


Format 

Description 

3 

VOCFS(I, J) 

4D15.8 

I = 1, 2, J = Local node number 
Nodal mean flow velocity values, I = 1 de- 
fines x-velocity and I = 2 defines r-velocity 

4 

ZVX(I) 

2D15.8 

I = Global node number 

Nodal mean flow x-component velocity val- 
ues 


ZVR(I) 

2D15.8 

I = Global node number 

Nodal mean flow r-component velocity val- 
ues 


n data descript i on of THE ACOUSTIC RAPIATIQN COD E 

1. Input data There are two different versions of the acoustic radiation pro- 
gram - PRATRADA and PRATRADB. PRATRADA is the version where 
the nodal acoustic pressure has been evaluated by an averaging technique while 
PRATRADB is the version where the acoustic pressure has been evaluated at the 
Gauss points and not at the nodal points. Therefore the postprocessors of the two 
versions are different but the finite element calculations are the same. The input 
data files for the acoustic radiation program PRATRADA and PRATRADB are 
files unit 5, 20 and 22. Data file 20 is the output of PRATMESH and contains 
the finite element mesh information. File 22 is the output of PRATV EL described 
before. The user input is in data file 5. The input data has been structured into 
cards. Each card of data begins on a new line and the input data is formatted. 
The acoustic radiation program is capable for running multiple cases. It contains 
the data for all the cases with an alphanumeric input separating the data for any 
two cases which determines whether the case is to be run or not. 
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Card 


Variable 


Format 


Description 


HDR(I) 


14A4 


I = 1, 14 

Character*4; it is the control for multiple 
cases; if HDR(l) = ‘stop’ the program stops 
execution and hence it is at the beginning of 
the data for each case 


NSYM 


110 


= 0, rectangular duct 
= 1, circular or annular duct 


ETAR 


F10.5 


Nondimensional frequency of the sound 
source (r? r = uR/c r ) 

u> = fan rotational speed in rad /sec 

R = reference duct radius at the fan face 
c r = freestream speed of sound outside the 
nacelle 


PRINTl 


15 


^ 0, Beginning row and column in stiffness 
matrices printed 

= 0, Beginning row and column in stiffness 
matrices not printed 


PRINT2 


15 


# 0, Final row and column in stiffness ma- 
trices printed 


= 0, Final row and column in stiffness ma- 
trices not printed 

NPLOT 15 > 0, level curves for the solution vector plot- 

ted 

= 0, plotting routine bypassed 
NCONT 15 Number of level curves to be plotted 
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Card 

3 


CMAXO 

F10.5 

Value of maximum level curve 

CMINO 

F10.5 

Value of minimum level curve 

ZAI(I) 

6F10.5 

1 = 1, NPOS; complex incident modal am- 
plitudes, real and imaginary parts 

NLINED 

15 

Number of acoustically lined elements on the 
inner surface of nacelle (if NLINED = 0, 
there is no lining and lining impedances are 
not required) 

MBEGIN 

15 

Element number (counted from fan face 
along the inner nacelle surface) on which the 
lining begins 

ZADM(I) 

6F10.5 

1=1, NLINED; admittances in the elements 
on the nacelle inner surface (complex val- 
ues), real and imaginary parts 


2, Output daU The output file for the acoustic radiation codes PRATRADA 
and PRAT RAD B is the file unit 6 which is well documented in itself and therefore 
is not described here. In addition to the saved files, seven unformatted scratch 
files - units 1, 2, 3, 4, 15, 16 and 17 are used. Files 15, 16, 17 are direct access 
and files 1, 2, 3, 4 are sequential access. Record lengths for the direct access files 
15, 16, and 17 are 36, 144 and 36 respectively. The number of records in each 
direct access file is the number of elements in the domain. In PRATRADB, the 
data for the acoustic pressure evaluated at the Gaussian quadrature points inside 
the elements has been written to file 6. The user may write it to a separate file 
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if he/she wishes to. It may then be used by some convenient plotting package 
like TEC PLOT to plot the sound pressure levels in the domain. The acoustic 
pressure data at the Gauss points is written down element by element. In each 
element there are four Gauss points. Therefore, the number of records for the 
acoustic pressure data is four times the total number of elements in the domain. 
The description of the acoustic pressure data evaluated at Gauss points follows. 


Card 

Variable 

Format 

Description 

1 

XG (I, J, K) 

1X,E14.7 

I = element number; J, K = 1, 2; (J, K) 
refers to a particular Gauss point in the 2 
x 2 grid; x coordinate of the Gauss point 


RG(I, J, K) 

2X,E14.7 

I = element number; J, K = 1, 2; (J, K) 
refers to a particular Gauss point in the 2 
x 2 grid; r coordinate of the Gauss point 


PGAU(I, J, K) 

2X.E14.7 

I = element number; J, K = 1, 2; (J, K) 
refers to a particular Gauss point in the 2 
x 2 grid; acoustic pressure at the Gauss 
point 


el data description of the cubic sptjne INTERPO LATION 

PROGRAM 

1. Input data This program is for generating the input data for nacelle and 
centerbody geometry. The nacelle and centerbody geometry is generated by a 
spline curve fit procedure using x and r coordinates of enough points to define 
the shape of the outer nacelle, the inner nacelle, and the center body. Using the 
spline information, the surfaces of the nacelle are discretized into line elements 
whose end points are defined. The center node of the elements is created from 
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the knowledge of the end points. The curve fit is a natural cubic spline. The 
cubic spline interpolation program prepares the input data for the mesh generation 
program PRATMESH. Details of the scheme have been described in section III. 
The input data is of free format. Each card of data begins on a new line. 


Card I Variable 


NELU 


Description 


Number of line elements describing the upper sur- 
face of the nacelle; the number of this parameter is 
one less than the number of points being input to 
represent the surface (see FRACTU(I)) 


NELC1 Number of line elements describing the centerbody 
between fan face and centerbody tip; the number of 
this parameter is one less than the number of points 
being input to represent the surface (see FRACl(I)) 


NELC2 Number of line elements describing the centerline 
between centerbody tip and intersection of highlight 
circle with z-axis; the number of this parameter is 
one less than the number of points being input to 
represent the surface (see FRAC2(I)) 


XBAF z-coordinate of the intersection of the baffle with 
the upper surface of the nacelle 


XTIP z-coordinate of the nacelle tip 


YTIP r-coordinate of the nacelle tip 


XFAN z-coordinate of the fan face 
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Card 

Variable 

Description 

2 

XCB 

z-coordinate of the centerbody tip 

3 

FRACTU(I) 

1=1, NELU + 1; Fractional distance of the end 
nodes of the line elements on the upper nacelle sur- 
face sequenced from baffle surface onwards (fraction 
based on the entire length of the upper nacelle sur- 
face) 

4 

FRACl(I) 

1=1, NELCl + 1; Fractional distance of the end 
nodes of the line elements on the centerbody se- 
quenced from fan face onwards (fraction based on 
the entire length of the centerbody) 

5 

FRAC2(I) 

1=1, NELC2 + 1; Fractional distance of the end 
nodes of the line elements on the centerline se- 
quenced from the centerbody tip to the highlight 
circle (fraction is based on this distance along the 
centerline) 


2. Output data The output of this program are the variables ADSHPU(I, J, K), 
ADSHPI(I, J, K), CBSHPI(I, J, K) (see input data description of MESHGEN ) 
written onto output file unit 7 with the same format as the corresponding data in 
input file 5 of MESHGEN. In a typical application, the spline program output 
file is imported into data file 5 for MESHGEN . A listing of the program follows. 


F. CUBIC SPLINE INTERPOLATION PROGRAM LISTING 

c******program for generating the input data for nacelle geometry***** 
implicit real*8 (a-h,o-z) 

dimension xnodu(200),yhodu(200),xnodi(200),ynodi(200),xnodc(200) 
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dimension ynodc(200) ,adshpu(100,3,2) ,adshpi( 100,3,2) 
dimension adshpc(100,3,2),fractu(200) 1 fracl(100),frac2(100) 
c xbaf = x coordinate of point where baffle starts 
c xtip = x coordinate of nacelle tip 

c ytip = y coordinate of nacelle tip 

c xfan = x coordinate of fan face 
c xcb = x coordinate of the end of centerbody 
read (5,*) nelu, nelcl, nelc2 
nelc = nelcl + nelc2 
nodeu = 2*nelu + 1 
nodec = 2* nelc + 1 
nodecl = 2*nelcl + 1 
nodec2 = 2*nelc2 + 1 
read(5,*) xbaf, xtip, ytip, xfan, xcb 
read(5,*) (fractu(i),i=l,nelu+l) 
read(5,*) (fracl(i),i=l,nelcl+l) 
read(5,*) (frac2(i),i=l,nelc2+l) 

c generating x coordinate nodal points on the upper nacelle surface 

do 20 i = 1, nelu 

xnodu(2*i-l) = xbaf + fractu(i)*(xtip - xbaf) 
xnodu(2*‘H-l) = xbaf + fractu(i+l)*(xtip - xbaf) 
xnodu(2*i) — (xnodu(2*i+l) 4- xnodu(2*i-l))/2.0d0 
20 continue 

call spline(nodeu,xnodu,ynodu) 

do 30 i = 1, nelu 
adshpu(i,l,l) = xnodu(2*i-l) 
adshpu(i,l,2) = ynodu(2*i-l) 
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adshpu(i,2,l) = xnodu(2*i) 
adshpu(i,2,2) = ynodu(2*i) 
adshpu(i,3,l) = xnodu(2*i+l) 
adshpu(i,3,2) = ynodu(2*i+l) 

write(7,1001) adshpu(i,l,l),adshpu(i,l,2),adshpu(i,2,l), 
&z adshpu(i,2,2),adshpu(i,3,l).adshpu(i,3,2) 

30 continue 

c generate the highlight circle 

c center of highlight circle 

chc = xtip - ytip 

c radius of highlight circle 

radhc = dsqrt(2.0d0*ytip*ytip) 
c. .intercept of highlight circle with the x-axis 
xhci = chc + radhc 

c generating x coordinate nodal points on the center body 

do 25 i = 1, nelcl 

xnodc(2*i-l) = xfan + fracl(i)*(xcb - xfan) 
xnodc(2*i+l) = xfan + fracl(i+l)*(xcb - xfan) 
xnodc(2*i) = (xnodc(2*i+l) + xnodc(2*i-l))/2.0d0 
write(6,*) ’xnodc(\2*i+l,’)=\ xnodc(2*i+l) 

25 continue 

c generating x coordinate nodal points on the centerline 

11 = 0 

do 31 i = nelcl 1, nelc 
11 = 11 + 1 

xnodc(2*i-l) = xcb + frac2(ll)*(xhci - xcb) 
xnodc(2*i+l) = xcb + frac2(ll+l)*(xhci - xcb) 
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xnodc(2*i) = (xnodc(2*i+l) -r xnodc(2*i-l))/2.0d0 
write(6,*) ’xnodc(’,2*i - 1,’)=’, xnodc(2*i- 1) 

31 continue 

call spline(nodec,xnodc,ynodc) 

do 35 i = 1, nelc 
adshpc(i,l,l) = xnodc(2*i-l) 
adshpc(i,l,2) = ynodc(2*i-l) 
adshpc(i,2,l) = xnodc(2*i) 
adshpc(i,2,2) = ynodc(2*i) 
adshpc(i,3,l) = xnodc(2*i+l) 
adshpc(i,3,2) = ynodc(2*i+l) 

35 continue 
c 

c.... compute a new set of fractions for the nodal points on the lower 

c... .surface of the nacelle which have the same fractions as the 

c.. corresponding nodes on the c.b & centerline based on the entire length 

dist = xtip - xfan 
do 75 k = 1, nodec 

fracti = (xnodc(k) - xfan)/(xhci - xfan) 
xnodi(k) = xfan + fracti* dist 
75 continue 

c do 76 k = no dec 1 + 1, nodec 
c fracti(k) = frac2(k)*(xhci * xcb)/(xhci - xfan) 
c xnodi(k) = xfan + fracti(k) dist 
c 76 continue 

call spline(nodec,xnodi,ynodi) 

do 85 i = 1, nelc 
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adshpi(i,l,l) = xnodi(2*i-l) 
adshpi(i,l)2) = ynodi(2*i-l) 
adshpi(i,2,l) = xnodi(2*i) 
adshpi(i,2,2) = ynodi(2*i) 
adshpi(i,3,l) = xnodi(2*i+l) 
adshpi(i,3,2) = ynodi(2*i+l) 

write(7,1001) adshpi(i,l,l),adshpi(i,l,2),adshpi(i,2,l), 

& adshpi(i,2,2),adshpi(i,3,l),adshpi(i,3,2) 

85 continue 

do 334 i = 1, nelc 

write(7,1001) adshpc(i,l,l)>adshpc(i,l,2),adshpc(i,2,l), 

& adshpc(i,2,2),adshpc(i,3,l),adshpc(i,3,2) 

334 continue 
1001 format (6fl0.4) 
stop 
end 

c ********************************************************************* 
c to fit a curve through a set of points using cubic spline 
c interpolation. 

c *************** ****************************************************** 
subroutine spline(nnode,pt,spl) 
implicit real*8 (a-h,o-z) 

dimension x(100) ,f( 100) ,b( 100) ,ed( 100) ,eu( 100) ,el( 100) ,dfp ( 100) 
dimension pt(200), spl(200) 

c read the data points from the data file 

read(5,*) n 

read(5,*) (x(i),f(i),i=l.n+l) 
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do 10 i = 2, n 
if(i.eq.2) go to 6 
el(i-l) = x(i) - x(i-l) 

6 ed(i-l) = 2.0d+00*(x(i+l) - x(i-l)) 
if(i.eq.n) go to 7 

eu(i-l) = x(i+l) - x(i) 

7 b(i-l) = 6.0d+00*(f(i+l) - f(i))/(x(i+l) - x(i)) 
k + 6.0d+00*(f(i-l) - f(i))/(x(i) - x(i-l)) 

10 continue 
nl = n - 1 

call tridag(nl,ed,eu,el,b) 
write(6,*) ’the solution is’ 
write(6,*) (b(i),i=l,n-l) 
dfp(l) = O.OdO 
dfp(n+l) = O.OdO 
do 70 i = 1, n-1 
dfp(i+l) = b(i) 

70 continue 

do 75 k = 1, nnode 
do 80 i — 1, n 

if(pt(k).gt.x(i).and.pt(k).lt.x(i+l))then 
spl(k) = dfp(i)*((x(i+l)-pt(k))**3)/(6.0d0*(x(i+l)-x(i))) 
+ dfp(i+l)*((pt(k)-x(i))**3)/(6.0d0*(x(i+l)-x(i))) 
k + (f(i)/(x(i+l)-x(i>) 

k - dfp(i)*(x(i+l)-x(i))/6.0d0)*(x(i+l)-pt(k)) 

k + (f(i+l)/(x(i+l)-x(i)) - dfp(i+l)*(x(i+l)-x(i))/6.0d0) 
k *(pt(k)-x(i)) 
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else if(pt(k).eq.x(i)) then 
spl(k) = f(i) 

else if(pt(k).eq.x(i+l)) then 

spl(k) = f(i+l) 

else if(pt(k).gt.x(n+l)) then 

spl(k) = O.OdO 

endif 

80 continue 
75 continue 

write(6,*) ’the calculated value is’ 
write(6,*) (spl(k),k=l,nnode) 
return 
end 
c 

subroutine tridag(nl,ed,eu,el,b) 

implicit real*8 (a-h,o-z) 

dimension ed(100) ,eu (100) ,el (100) ,b (100) 

m = nl - 1 

do 1040 i = l,m 

fa = el(i+l)/ed(i) 

ed(i+l) = ed(i+l) - fa*eu(i) 

b(i+l) = b(i+l) - fa*b(i) 

1040 continue 

b(nl) = b(nl)/ed(nl) 
do 1070 i = l,m 

b(nl-i) = (b(nl-i) - eu(nl-i)*b(nl-i+l))/ed(nl-i) 
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1070 


continue 

return 

end 
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APPENDIX B 


MATHEMATICAL NOTATIONS 


rp-.n — ► R 

xl > : 0 — ♦ R 2 
y\) : n — ► C 2 


|0,r o 


3 

V 

u 

H' 


= ^ is a real valued function defined on the domain H which 
is one-dimensional 

= -0 is a. real valued function defined on the domain H which 
is two-dimensional 

= xJj is a complex valued function defined on the domain Cl which 
is two-dimensional 

= an interval between and including the values of 0 and r a 
= such that 
= for every 

= symbol representing union of two sets 
= Hilbert space; the space of functions which are square 
integrable, which, in other words, implies that the 
functions are continuous, but their derivatives 
are piecewise continuous 
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ACOUSTIC RADIATION CODE 
INSTALLATION NOTES 


INTRODUCTORY REMARKS 

The Acoustic Radiation Code is a FORTRAN program used to study far field radiation from 
turbofan engines. It was developed for the IBM (tm) mainframe at the University of Missouri- 
Rolla, under the direction of Professor Walter Eversman. It has subsequently been modified at 
Hamilton Standard to run on Sun (tm) and Silicon Graphics Iris (tm) UNIX (tm) workstations. 

The program consists of five separate modules. These are run, one after the other, in the order 
of their listing in the next section. However, PRATPREH, the first of these, may not always 
be used. 

The purpose of these notes is to assist users in the installation of the code on either of the two 
above-mentioned workstations. In the pages that follow, there is a brief description of the 
modules making up the program and then brief descriptions of how to compile, run, and test 
these modules. 

For further details, refer to Appendix A. 


ACOUSTIC RADIATION CODE MQDIITES 
Five modules have been provided for installation: 

PRATPREH 

► Generates cards 3-5 of the PRATMESH.INP input to PRATMESHH (i.e., 
the coordinates of the element nodes for the upper nacelle, center body and 
lower nacelle) 
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PRATMESHH 


► Generates the finite element mesh used for both flow and acoustic solutions 

► Calculates the duct eigenvalues and eigenfunctions used for acoustic 
calculations 


PRATFLOWH 

► Obtains the potential flow solutions 


PRATVELH 


* Provides a superposition of the solutions above with a uniform mean flow 
to give the final flow needed by PRATRADH 

PRATRADH 

► Generates the acoustic solution for two-dimensional or cylindrically 
symmetric nacelles 

Note that the equivalent University of Missouri-Rolla IBM (tm) mainframe versions lof these 
modules are designated PRATPRE, PRATMESH, PRATFLOW, PRATVEL, and PRATRADA. 
An "H" for Hamilton has been added to the names of the workstation versions (or in the case 
of PRATRADA, the final "A" has been changed to "H"), to distinguish between the two 

versions. 


TO COMPILE 

PRATPREH 

f77 -o pratpreh -O pratpreh.f ( Sun ) 

f77 -o pratpreh -O -old_rl pratpreh.f (SGI) 

PRATMESHH 

ni -o pratmeshh -O pratmeshh.f (Sun) 

ni -o pratmeshh -0 -old_rl pratmeshh.f (SGI) 
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PRATFLOWH 


f77 -o pratflowh -O pratflowh.f (Sun) 

f77 -o pratflowh -O -old_rl pratflowh.f (SGI) 

PRATVELH 


f77 -o pratvelh 
f77 -o pratvelh 
PRATRADH 

f77 -o pratradh 
f77 -o pratradh 


-O pratvelh. f (S un ) 

-O -old_rl pratvelh. f (SGI) 

-O pratradh. f (Sun) 

-O -Olimit 1100 -old_rl pratradh.f (SGI) 


TO RUN 

PRATPREH 


(pratpreh < pratpre.inp > pratpre.out) >& pratpre.err & 

(fort.5) (fort. 6) (stderr) 

Files Generated: 


* fort. 7 - Cards 3-5 for pratmesh.inp 

* pratpre.out - Numeric output data 

* pratpre.err - System error messages 

PRATMESHH 


(3.3 Kb) 
(9 Kb) 


(pratmeshh < pratmesh.inp > pratmesh.out) >& pratmesh.err & 

(fort.5) (fort. 6) (stderr) 


Files Generated: 


* fort. 14 

* fort. 20 

* pratmesh.out 

* pratmesh.err 


- PostScript plot file (450 Kb) 

- Needed as input by all other modules (732 Kb) 

- Numeric output data (360 Kb) 

- System error messages 
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PRATFLOWH 


(pratflowh < pratflow.inp > pratflow.out) >& pratflow.err & 
(fort.5) (fort.6) (stderr) 

Note: fort.20, generated by PRATMESHH, must be 

available as UNIT =20 input. 


Files Generated: 


* fort. 2-4, 8 

& 5 tmp files 

* fort. 14 

* fort. 21 

* pratflow.out 

* pratflow.err 


- Work files; delete 
the fort. 2-4, 8 ones 
at end of run 

- PostScript plot file 

- Needed as input by PRATVELH 

- Numeric output data 

- System error messages 


(6.5 Mb total) 


(1.6 Mb) 
(195 Kb) 
(196 Kb) 


PRATVELH 

(pratvelh < pratvel.inp > pratvel.out) >& pratvel.err & 

(fort.5) (fort.6) (stderr) 

Note: fort.20, generated by PRATMESHH, and fort.21, 

generated by PRATFLOWH, must be available, 
respectively, as UNIT =20 and UNIT =21 input. 


Files Generated: 


* 4 tmp files 

* fort. 14 

* fort. 22 

* pratvel.out 

* pratvel.err 


- Work files 

- PostScript plot file 

- Needed as input by PRATRADH 

- Numeric output data 

- System error messages 


(443 Kb total) 
(513 Kb) 

(700 Kb) 

(13 Kb) 


PRATRADH 

(pratradh < pratrad.inp > pratrad.out) >& pratrad.err & 

(fort.5) (fort.6) (stderr) 


Note: fort.20, generated by PRATMESHH, and fort. 22, 

generated by PRATVELH, must be available, 
respectively, as UN1T=20 and UN1T=22 input. 
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Files Generated: 


* fort. 3, 4 & 

5 tmp files 

* fort. 14 

* fort.23 (and 
fort. 24-27, 

one for each 
additional case 
that is run) 

* pratrad.out 

* pra trad, err 


Work files; delete the 
fort 3, 4 ones at end of run 
PostScript plot file 

Used previously for additional plotting; 
not needed for the present 
setup, so can be deleted 


Numeric output data 
System error messages 


(70 Mb total) 

(1 Mb per case) 
(15 Kb per case) 


(107 Kb per case) 


Note that the plot files, fort. 14, are easily plotted using a standard PostScript printer. They also 
can be previewed on a monitor if a PostScript viewing utility is available. 
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TEST CASE 


A test case has been provided to check each of the five Acoustic Radiation Code modules. For 
test purposes PRATPREH can be run independently. However, the other programs should be 
run in sequence with the output (i.e., fort. 20, fort.21, fort. 22) of one feeding into the next. The 
test input Files to use have the same names as those used in the "To Run" instructions above 
(i.e., pratpre.inp, pratmesh.inp, pratflow.inp, pratvel.inp, pratrad.inp). 

Output data (fort. 7) from the PRATPRE test run should be as follows: 


0.0000 
0.0450 
0.0900 
0.1350 
0.1800 
0.2250 
0.2700 
0.3150 
0.3600 
0.4050 
0.4500 
0.4950 
0.5400 
0 .5850 
0.6300 
0.6750 
0.7200 
0.7650 
0.8100 
0.8550 
0.8775 
0 .0000 
0 .0415 
0.0830 
0.1245 
0 . 1660 
0.2075 
0.2490 
0.2905 
0.3319 
0 .3717 
0 .4115 
0.4512 
0.4910 
0.5308 
0.5705 
0 .6103 
0 . 6 S 01 
0.6898 
0.7298 
0.7693 
0.8091 
0.8489 
0 . 8688 
0.0000 
0.0625 
0. 12S0 

0 . 1875 
0.2500 
0.3125 
0.3750 
0 .4375 
0.5000 
0.5599 
0.6196 
0 . 6797 
0.7396 
0.7995 
0 . 8594 
0.9193 
0.9792 
1.0390 

1 . 0989 
1 . 1588 
1.2187 
1 .2786 
1.3086 


1.2000 
1.2000 
1.2000 
1.2000 
1.2000 
1.2000 
1.2000 
1.2000 
1.2000 
1.2000 
1.2000 
1.2000 
1.2000 
1.2000 
1.1995 
1.1968 
1.1917 
1.1835 
1.1714 
1.1527 
1.1381 
1.0000 
1.0000 
1.0000 
1.0000 
1.0000 
1.0000 
1.0000 
1.0000 
1.0000 
1.0000 
1.0000 
1.0000 
1.0000 
1.0000 
1.0000 
1.0001 
1.0014 
1.0046 
1.0098 
1.0175 
1.0283 
1.0441 
1.0556 
0.3000 
0.2818 
0.2601 
0.2328 
0.2003 
0.1619 
. 0.1166 
0.0635 
0.0000 
0.0000 
0.0000 
0.0000 
0 .0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 


0.0225 

0.0675 

0.1125 

0.1575 

0.2025 

0.2475 

0.2925 

0.3375 

0.3825 

0 . 4275 

0.4725 

0.5175 

0.5625 

0.6075 

0 . 6 S 25 

0.6975 

0 . 742 S 

0.7875 

0.8325 

0.8662 

0 . B 887 

0.0207 

0.0622 

0.1037 

0.1452 

0 . 1867 

0.2282 

0.2697 

0.3112 

0.3518 

0.3916 

0.4314 

0.4711 

0.5109 

0.5506 

0.5904 

0.6302 

0.6699 

0.7097 

0.7495 

0.7892 

0.8290 

0.8588 

0.8844 

0.0312 

0.0938 

0 . 1562 

0.2188 

0.2812 

0.3438 

0.4062 

0.4688 

0.5299 

0 .5898 

0.6497 

0.7096 

0.7695 

0.8294 

0.8893 

0.9492 

1.0091 

1.0690 

1 . 1289 

1.1888 

1.2487 

1.2936 

1.3321 


1.2000 
1.2000 
1.2000 
1.2000 
1.2000 
1.2000 
1.2000 
1.2000 
1.2000 
1.2000 
1.2000 
1.2000 
1.2000 
1.1999 
1.1985 
1.1946 
1.1880 
1.1781 
1.1632 
1.1460 
1.1269 
1.0000 
1.0000 
1.0000 
1.0000 
1.0000 
1.0000 
1.0000 
1.0000 
1.0000 
1.0000 
1.0000 
1.0000 
1.0000 
1.0000 
1.0000 
1.0006 
1.0027 
1.0069 
1.0133 
1.0224 
1.0354 
1 . 0495 
1.0685 
0.2911 
0.2716 
0.2471 
0.2172 
0.1819 
0.1402 
0.0911 
0.0332 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0 .0000 
0.0000 
0.0000 
0.0000 
0.0000 


0.0450 
0.0900 
0.1350 
0.1800 
0.2250 
0.2700 
0.3150 
0.3600 
0.4050 
0.4500 
0.4950 
0.5400 
0.5850 
0.6300 
0.6750 
0.7200 
0.7550 
0.8100 
0.8550 
0.8775 
0.9000 
0.0415 
0.0830 
0.1245 
0.1660 
0 .2075 
0.2490 
0 . 2905 
0.3319 
0.3717 
0.4115 
0.4512 
0.4910 
0.5308 
0.5705 
0.6103 
0.6501 
0.6898 
0.7296 
0.7693 
0.8091 
0.8489 
0.8688 
0.9000 
0.0625 
0.1250 
0.1875 
0.2500 
0.3125 
0.3750 
0.4375 
0.5000 
0.5599 
0.6198 
0.6797 
0 .7396 
0.7995 
0.8594 
0.9193 
0.9792 
1.0390 
1.0989 
1 . 1588 
1.2187 
1.2786 
1 .3086 
1.3556 


1.2000 
1.2000 
1.2000 
1.2000 
1.2000 
1.2000 
1.2000 
1.2000 
1.2000 
1.2000 
1.2000 
1.2000 
1.2000 
1.1995 
1.1968 
1.1917 
1.1835 
1.1714 
1.1527 
1.1381 
1.1000 
1.0000 
1.0000 
1.0000 
1.0000 
1.0000 
1.0000 
1.0000 
1.0000 
1.0000 
l .0000 
1.0000 
1.0000 
1.0000 
1.0000 
1.0001 
1.0014 
1.0046 
1.0098 
1.0175 
1.0283 
1.0441 
1.0556 
1 .1000 
0.2818 
0.2601 
0.2328 
0.2003 
0.1619 
0 . 1166 
0.0635 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
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Output plots from runs of the remaining modules should match those shown below. Note that 
plots here are reduced in size. The actual ones will be 8V2" x 11". 


PRATMESHH 
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